ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/remailleur.cpp
Revision: 147
Committed: Wed Aug 27 20:35:33 2008 UTC (16 years, 8 months ago) by francois
Original Path: magic/lib/mailleur/mailleur/src/remailleur.cpp
File size: 44184 byte(s)
Log Message:
Remailleur - reprise du doc dans l environnement MAGiC

File Contents

# User Rev Content
1 francois 147 //---------------------------------------------------------------------------
2    
3     #include "gestionversion.h"
4     #include <math.h>
5     #include "remailleur.h"
6     #include "mg_gestionnaire.h"
7     #include "mailleur3d.h"
8     #include "mailleur0d.h"
9     #include "mailleur1d.h"
10     #include "mailleur2d.h"
11     #include "ot_mathematique.h"
12     #include "mg_segment_frontiere.h"
13     //---------------------------------------------------------------------------
14     #pragma package(smart_init)
15    
16    
17    
18    
19    
20     REMAILLEUR::REMAILLEUR(class MG_GESTIONNAIRE* g1,MG_GESTIONNAIRE* g2,MG_GEOMETRIE* geo1,MG_GEOMETRIE* geo2,class FCT_TAILLE* fct_taille1,FCT_TAILLE* fct_taille2,MG_MAILLAGE* mori,MG_MAILLAGE* mmodi,VCT_COMPARAISON_RESULTAT& cmp):gestorigine(g1),gestmodifie(g2),carteori(fct_taille1),cartemod(fct_taille2),rescmp(cmp),maiorigine(mori),maimodifie(mmodi),geoorigine(geo1),geomodifie(geo2)
21     {
22     int nbface1=geomodifie->get_nb_mg_face();
23     quad=new TPL_QUADTREE<MG_SEGMENT_FRONTIERE*,MG_NOEUD*>*[nbface1];
24     int nbvolume=geomodifie->get_nb_mg_volume();
25     lsttrifront=new TPL_LISTE_ENTITE<MG_TRIANGLE*>[nbvolume];
26     }
27    
28     REMAILLEUR::~REMAILLEUR()
29     {
30     int nbface1=geomodifie->get_nb_mg_face();
31     for (int i=0;i<nbface1;i++) delete quad[i];
32     delete [] quad;
33     TPL_MAP_ENTITE<class MG_SEGMENT_FRONTIERE*>::ITERATEUR it;
34     for (MG_SEGMENT_FRONTIERE* seg=lstsegfront.get_premier(it);seg!=NULL;seg=lstsegfront.get_suivant(it) ) delete seg;
35     delete octree;
36     delete [] lsttrifront;
37     }
38    
39    
40     void REMAILLEUR::maille(void)
41     {
42     maille(20);
43     }
44    
45     void REMAILLEUR::maille(int etape)
46     {
47     //etape 0 initialisation
48     initialise_octree();
49     //etape 1 Destruction des entites autour des entites disparues
50     if (etape<1) return;
51     affiche(" Destruction autour des disparitions");
52     TPL_MAP_ENTITE<MG_NOEUD*> lstnoeuddetruire;
53     int nb=rescmp.get_nb_liste_topologie(ORIGINE_DISPARUE);
54     for (int i=0;i<nb;i++)
55     {
56     MG_ELEMENT_TOPOLOGIQUE* ele=rescmp.get_liste_topologie(ORIGINE_DISPARUE,i);
57     int dim=ele->get_dimension();
58     if (dim==0)
59     {
60     MG_SOMMET* som=(MG_SOMMET*)ele;
61     int nblien=som->get_lien_maillage()->get_nb();
62     for (int j=0;j<nblien;j++)
63     {
64     MG_NOEUD* no=(MG_NOEUD*)som->get_lien_maillage()->get(j);
65     if (maiorigine->get_mg_noeudid(no->get_id())!=NULL) lstnoeuddetruire.ajouter(no);
66     }
67     }
68     if (dim==1)
69     {
70     MG_ARETE* are=(MG_ARETE*)ele;
71     int nblien=are->get_lien_maillage()->get_nb();
72     for (int j=0;j<nblien;j++)
73     {
74     MG_SEGMENT* seg=(MG_SEGMENT*)are->get_lien_maillage()->get(j);
75     if (maiorigine->get_mg_segmentid(seg->get_id())!=NULL)
76     {
77     lstnoeuddetruire.ajouter(seg->get_noeud1());
78     lstnoeuddetruire.ajouter(seg->get_noeud2());
79     }
80     }
81     }
82     if (dim==2)
83     {
84     MG_FACE* face=(MG_FACE*)ele;
85     int nblien=face->get_lien_maillage()->get_nb();
86     for (int j=0;j<nblien;j++)
87     {
88     MG_TRIANGLE* tri=(MG_TRIANGLE*)face->get_lien_maillage()->get(j);
89     if (maiorigine->get_mg_triangleid(tri->get_id())!=NULL)
90     {
91     lstnoeuddetruire.ajouter(tri->get_noeud1());
92     lstnoeuddetruire.ajouter(tri->get_noeud2());
93     lstnoeuddetruire.ajouter(tri->get_noeud3());
94     }
95     }
96     }
97     }
98     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR it;
99     for (MG_NOEUD* no=lstnoeuddetruire.get_premier(it);no!=NULL;no=lstnoeuddetruire.get_suivant(it))
100     {
101     double ereelle=0.;
102     int nbseg=no->get_lien_segment()->get_nb();
103     for (int i=0;i<nbseg;i++)
104     ereelle=ereelle+no->get_lien_segment()->get(i)->get_longueur();
105     ereelle=ereelle/nbseg;
106     detruit_noeud(no,2*ereelle);
107     }
108     //etape 2 : transfert des noeuds de sommets;
109     if (etape<2) return;
110     affiche(" Transfert sommet");
111     nb=rescmp.get_nb_liste_topologie(ORIGINE_CONSERVEE);
112     for (int i=0;i<nb;i++)
113     {
114     MG_ELEMENT_TOPOLOGIQUE* ele=rescmp.get_liste_topologie(ORIGINE_CONSERVEE,i);
115     int dim=ele->get_dimension();
116     if (dim==0)
117     {
118     MG_SOMMET* som=(MG_SOMMET*)ele;
119     int nblien=som->get_lien_maillage()->get_nb();
120     for (int j=0;j<nblien;j++)
121     {
122     MG_NOEUD* no=(MG_NOEUD*)som->get_lien_maillage()->get(j);
123     if ((maiorigine->get_mg_noeudid(no->get_id())!=NULL) && (no->get_nouveau_numero()==CONSERVE))
124     {
125     MG_SOMMET* nvsom=geomodifie->get_mg_sommetid(rescmp.get_liste_correspondance(som->get_id()));
126     transfert_noeud(no,nvsom);
127     }
128     }
129    
130     }
131     }
132     // etape 3 maillage 0d
133     if (etape<3) return;
134     affiche(" MAILLAGE 0D");
135     MAILLEUR0D m0d(maimodifie,geomodifie);
136     m0d.adapte();
137     // etape 4 destruction autour des nouveaux sommets
138     if (etape<4) return;
139     affiche(" Destruction autour des nouveaux sommets");
140     lstnoeuddetruire.vide();
141     nb=rescmp.get_nb_liste_topologie(MODIFIE_APPARUE);
142     for (int i=0;i<nb;i++)
143     {
144     MG_ELEMENT_TOPOLOGIQUE* ele=rescmp.get_liste_topologie(MODIFIE_APPARUE,i);
145     int dim=ele->get_dimension();
146     if (dim==0)
147     {
148     MG_SOMMET* som=(MG_SOMMET*)ele;
149     int nblien=som->get_lien_maillage()->get_nb();
150     for (int j=0;j<nblien;j++)
151     {
152     MG_NOEUD* no=(MG_NOEUD*)som->get_lien_maillage()->get(j);
153     if (maimodifie->get_mg_noeudid(no->get_id())!=NULL) lstnoeuddetruire.ajouter(no);
154     }
155     }
156     }
157     for (MG_NOEUD* no=lstnoeuddetruire.get_premier(it);no!=NULL;no=lstnoeuddetruire.get_suivant(it))
158     {
159     detruit_noeud(no,-1,1);
160     }
161     // etape 5 transferts des segments d arete
162     if (etape<5) return;
163     /*affiche(" Transfert arete");
164     nb=rescmp.get_nb_liste_topologie(ORIGINE_CONSERVEE);
165     for (int i=0;i<nb;i++)
166     {
167     MG_ELEMENT_TOPOLOGIQUE* ele=rescmp.get_liste_topologie(ORIGINE_CONSERVEE,i);
168     int dim=ele->get_dimension();
169     if (dim==1)
170     {
171     MG_ARETE* are=(MG_ARETE*)ele;
172     int nblien=are->get_lien_maillage()->get_nb();
173     for (int j=0;j<nblien;j++)
174     {
175     MG_SEGMENT* seg=(MG_SEGMENT*)are->get_lien_maillage()->get(j);
176     if ((maiorigine->get_mg_segmentid(seg->get_id())!=NULL))
177     {
178     MG_NOEUD* n1=seg->get_noeud1();
179     MG_NOEUD* n2=seg->get_noeud2();
180     if (n1->get_nouveau_numero()==CONSERVE)
181     if (n2->get_nouveau_numero()==CONSERVE)
182     {
183     MG_ARETE* nvare=geomodifie->get_mg_areteid(rescmp.get_liste_correspondance(are->get_id()));
184     MG_NOEUD* nv1=get_liste_correspondance(n1);
185     MG_NOEUD* nv2=get_liste_correspondance(n2);
186     if (nv1==NULL) nv1=transfert_noeud(n1,nvare);
187     if (nv2==NULL) nv2=transfert_noeud(n2,nvare);
188     MG_SEGMENT *seg=new MG_SEGMENT(nvare,nv1,nv2,IMPOSE);
189     maimodifie->ajouter_mg_segment(seg);
190     }
191     }
192     }
193    
194     }
195     }
196     */
197     affiche(" Creation des blocs de maille 1D");
198     int nbarete=geomodifie->get_nb_mg_arete();
199     nb=rescmp.get_nb_liste_topologie(ORIGINE_CONSERVEE);
200     for (int i=0;i<nb;i++)
201     {
202     MG_ELEMENT_TOPOLOGIQUE* ele=rescmp.get_liste_topologie(ORIGINE_CONSERVEE,i);
203     int dim=ele->get_dimension();
204     if (dim==1)
205     {
206     multimap<unsigned long,MG_SEGMENT*,less <unsigned long> > lstseg;
207     MG_ARETE* arete=(MG_ARETE*)ele;
208     TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR ittet;
209     for (MG_ELEMENT_MAILLAGE* ele=arete->get_lien_maillage()->get_premier(ittet);ele!=NULL;ele=arete->get_lien_maillage()->get_suivant(ittet))
210     {
211     MG_SEGMENT* seg=(MG_SEGMENT*)ele;
212     if ((maiorigine->get_mg_segmentid(seg->get_id())!=NULL))
213     {
214     MG_NOEUD* n1=seg->get_noeud1();
215     MG_NOEUD* n2=seg->get_noeud2();
216     if (n1->get_nouveau_numero()==CONSERVE)
217     if (n2->get_nouveau_numero()==CONSERVE)
218     {
219     seg->change_nouveau_numero(NONTRAITE);
220     pair<unsigned long,MG_SEGMENT*> tmp(n1->get_id(),seg);
221     lstseg.insert(tmp);
222     pair<unsigned long,MG_SEGMENT*> tmp2(n2->get_id(),seg);
223     lstseg.insert(tmp2);
224     }
225     }
226     }
227    
228     cree_bloc_maille(lstseg,arete);
229     }
230     }
231     // etape 6 Maillage des aretes
232     if (etape<6) return;
233     affiche(" MAILLAGE 1D");
234     int nbare=geomodifie->get_nb_mg_arete();
235     for (int i=0;i<nbare;i++)
236     {
237     MG_ARETE* arete=geomodifie->get_mg_arete(i);
238     TPL_LISTE_ENTITE<double> param;
239     int typeface=arete->get_courbe()->get_type_geometrique(param);
240     multimap<unsigned long,BLOC_MAILLE_1D*,less <unsigned long> >::iterator it=lstb1d.begin();
241     while (it!=lstb1d.end())
242     {
243     BLOC_MAILLE_1D *b1d=(*it).second;
244     if ((b1d->etat==NONATTACHE) && (b1d->type==typeface))
245     {
246     MG_NOEUD* no1=b1d->lst[0]->get_noeud1();
247     MG_NOEUD* no2=b1d->lst[0]->get_noeud2();
248     double xyz1[3],xyz2[3],xyz[3];
249     transfert_coord(no1->get_coord(),xyz1);
250     transfert_coord(no2->get_coord(),xyz2);
251     xyz[0]=0.5*(xyz1[0]+xyz2[0]);
252     xyz[1]=0.5*(xyz1[1]+xyz2[1]);
253     xyz[2]=0.5*(xyz1[2]+xyz2[2]);
254     double t;
255     arete->inverser(t,xyz);
256     double xyztmp[3];
257     arete->evaluer(t,xyztmp);
258     OT_VECTEUR_3D vec(xyz,xyztmp);
259     double metrique[9];
260     cartemod->evaluer(xyz,metrique);
261     double eps=(1./sqrt(metrique[0])*1e-6);
262     if (vec.get_longueur()<eps)
263     {
264     if (arete->get_courbe()->est_periodique()==1)
265     if (t<arete->get_tmin()) t=t+arete->get_courbe()->get_periode();
266     if ((t>=arete->get_tmin()) && (t<=arete->get_tmax()))
267     {
268     int nbseg=b1d->lst.size();
269     b1d->etat=ATTACHE;
270     for (int i=0;i<nbseg;i++)
271     {
272     MG_SEGMENT *seg=b1d->lst[i];
273     MG_NOEUD* no1=seg->get_noeud1();
274     MG_NOEUD* no2=seg->get_noeud2();
275     double xyz1[3],xy2[3];
276     transfert_coord(no1->get_coord(),xyz1);
277     transfert_coord(no2->get_coord(),xyz2);
278     MG_NOEUD* nv1=get_liste_correspondance(no1);
279     MG_NOEUD* nv2=get_liste_correspondance(no2);
280     if (nv1==NULL) nv1=transfert_noeud(no1,arete);
281     if (nv2==NULL) nv2=transfert_noeud(no2,arete);
282     MG_SEGMENT *nvseg=new MG_SEGMENT(arete,nv1,nv2,IMPOSE);
283     maimodifie->ajouter_mg_segment(nvseg);
284     }
285     }
286     }
287     }
288     it++;
289     }
290     MAILLEUR1D m1d(maimodifie,geomodifie,arete,cartemod);
291     m1d.adapte();
292     }
293     cree_quadtree(maimodifie,quad,&lstsegfront);
294     // etape 7 Destruction autour des nouveaux segments
295     if (etape<7) return;
296     affiche(" Destruction autour des nouvelles aretes");
297     lstnoeuddetruire.vide();
298     nb=rescmp.get_nb_liste_topologie(MODIFIE_APPARUE);
299     for (int i=0;i<nb;i++)
300     {
301     MG_ELEMENT_TOPOLOGIQUE* ele=rescmp.get_liste_topologie(MODIFIE_APPARUE,i);
302     int dim=ele->get_dimension();
303     if (dim==1)
304     {
305     MG_ARETE* are=(MG_ARETE*)ele;
306     int nblien=are->get_lien_maillage()->get_nb();
307     for (int j=0;j<nblien;j++)
308     {
309     MG_SEGMENT* seg=(MG_SEGMENT*)are->get_lien_maillage()->get(j);
310     if (maimodifie->get_mg_segmentid(seg->get_id())!=NULL)
311     {
312     lstnoeuddetruire.ajouter(seg->get_noeud1());
313     lstnoeuddetruire.ajouter(seg->get_noeud2());
314     }
315     }
316     }
317     }
318     for (MG_NOEUD* no=lstnoeuddetruire.get_premier(it);no!=NULL;no=lstnoeuddetruire.get_suivant(it))
319     {
320     detruit_noeud(no,-1,1);
321     }
322    
323     // etape 8 Creation des blocs de mailles 2D
324     if (etape<8) return;
325     affiche(" Creation des blocs de maille 2D");
326     int nbface=geomodifie->get_nb_mg_face();
327     nb=rescmp.get_nb_liste_topologie(ORIGINE_CONSERVEE);
328     for (int i=0;i<nb;i++)
329     {
330     MG_ELEMENT_TOPOLOGIQUE* ele=rescmp.get_liste_topologie(ORIGINE_CONSERVEE,i);
331     int dim=ele->get_dimension();
332     if (dim==2)
333     {
334     multimap<unsigned long,MG_TRIANGLE*,less <unsigned long> > lsttri;
335     MG_FACE* face=(MG_FACE*)ele;
336     TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR ittet;
337     for (MG_ELEMENT_MAILLAGE* ele=face->get_lien_maillage()->get_premier(ittet);ele!=NULL;ele=face->get_lien_maillage()->get_suivant(ittet))
338     {
339     MG_TRIANGLE* tri=(MG_TRIANGLE*)ele;
340     if ((maiorigine->get_mg_triangleid(tri->get_id())!=NULL))
341     {
342     MG_NOEUD* n1=tri->get_noeud1();
343     MG_NOEUD* n2=tri->get_noeud2();
344     MG_NOEUD* n3=tri->get_noeud3();
345     if (n1->get_nouveau_numero()==CONSERVE)
346     if (n2->get_nouveau_numero()==CONSERVE)
347     if (n3->get_nouveau_numero()==CONSERVE)
348     {
349     tri->change_nouveau_numero(NONTRAITE);
350     pair<unsigned long,MG_TRIANGLE*> tmp(n1->get_id(),tri);
351     lsttri.insert(tmp);
352     pair<unsigned long,MG_TRIANGLE*> tmp2(n2->get_id(),tri);
353     lsttri.insert(tmp2);
354     pair<unsigned long,MG_TRIANGLE*> tmp3(n3->get_id(),tri);
355     lsttri.insert(tmp3);
356     }
357     }
358     }
359    
360     cree_bloc_maille(lsttri,face);
361     }
362     }
363     // etape 9 MAILLAGE 2D
364     if (etape<9) return;
365     affiche(" MAILLAGE 2D");
366     nbface=geomodifie->get_nb_mg_face();
367     for (int i=0;i<nbface;i++)
368     {
369     TPL_LISTE_ENTITE<MG_TRIANGLE*> lsttricontraint;
370     MG_FACE* face=geomodifie->get_mg_face(i);
371     TPL_LISTE_ENTITE<double> param;
372     int typeface=face->get_surface()->get_type_geometrique(param);
373     multimap<unsigned long,BLOC_MAILLE_2D*,less <unsigned long> >::iterator it=lstb2d.begin();
374     while (it!=lstb2d.end())
375     {
376     BLOC_MAILLE_2D *b2d=(*it).second;
377     if ((b2d->etat==NONATTACHE) && (b2d->type==typeface))
378     {
379     MG_NOEUD* no1=b2d->lst[0]->get_noeud1();
380     MG_NOEUD* no2=b2d->lst[0]->get_noeud2();
381     MG_NOEUD* no3=b2d->lst[0]->get_noeud3();
382     double uv[2],uv1[2],uv2[2],uv3[2];
383     double xyz1[3],xyz2[3],xyz3[3];
384     transfert_coord(no1->get_coord(),xyz1);
385     transfert_coord(no2->get_coord(),xyz2);
386     transfert_coord(no3->get_coord(),xyz3);
387     face->inverser(uv1,xyz1);
388     face->inverser(uv2,xyz2);
389     face->inverser(uv3,xyz3);
390     double xyz1tmp[3];
391     face->evaluer(uv1,xyz1tmp);
392     OT_VECTEUR_3D vec(xyz1,xyz1tmp);
393     double metrique[9];
394     cartemod->evaluer(xyz1,metrique);
395     double eps=(1./sqrt(metrique[0])*1e-6);
396     if (vec.get_longueur()<eps)
397     {
398     OT_DECALAGE_PARAMETRE deca(face->get_surface()->get_periode_u(),face->get_surface()->get_periode_v());
399     double du=deca.calcul_decalage_parametre_u(uv1[0]);
400     double dv=deca.calcul_decalage_parametre_v(uv1[1]);
401     uv[0]=0.333333333333333333*(deca.decalage_parametre_u(uv1[0],du)+deca.decalage_parametre_u(uv2[0],du)+deca.decalage_parametre_u(uv3[0],du));
402     uv[1]=0.333333333333333333*(deca.decalage_parametre_v(uv1[1],dv)+deca.decalage_parametre_v(uv2[1],dv)+deca.decalage_parametre_v(uv3[1],dv));
403     uv[0]=deca.decalage_parametre_u(uv[0],-du);
404     uv[1]=deca.decalage_parametre_v(uv[1],-dv);
405     double dist=calcule_distance_contour_face(uv,face,i,quad);
406     if (dist>0.)
407     {
408     int nbtri=b2d->lst.size();
409     b2d->etat=ATTACHE;
410     for (int i=0;i<nbtri;i++)
411     {
412     MG_TRIANGLE *tri=b2d->lst[i];
413     MG_NOEUD* no1=tri->get_noeud1();
414     MG_NOEUD* no2=tri->get_noeud2();
415     MG_NOEUD* no3=tri->get_noeud3();
416     double xyz1[3],xy2[3],xyz3[3];
417     transfert_coord(no1->get_coord(),xyz1);
418     transfert_coord(no2->get_coord(),xyz2);
419     transfert_coord(no3->get_coord(),xyz3);
420     MG_NOEUD* nvno1=new MG_NOEUD(NULL,xyz1[0],xyz1[1],xyz1[2],IMPOSE);
421     MG_NOEUD* nvno2=new MG_NOEUD(NULL,xyz2[0],xyz2[1],xyz2[2],IMPOSE);
422     MG_NOEUD* nvno3=new MG_NOEUD(NULL,xyz3[0],xyz3[1],xyz3[2],IMPOSE);
423     MG_TRIANGLE* nvtri=new MG_TRIANGLE(NULL,nvno1,nvno2,nvno3,NULL,NULL,NULL,IMPOSE);
424     lsttricontraint.ajouter(nvtri);
425     }
426     }
427     }
428     }
429     it++;
430     }
431     MAILLEUR2D m2d(maimodifie,geomodifie,face,cartemod);
432     char mess[50];
433     sprintf(mess," Face %i identificateur %lu",i,face->get_id());
434     affiche(mess);
435     m2d.maille(face,NULL,&(lsttricontraint));
436     int nbtricon=lsttricontraint.get_nb();
437     for (int i=0;i<nbtricon;i++)
438     {
439     MG_TRIANGLE* tri=lsttricontraint.get(0);
440     MG_NOEUD* noeud1=tri->get_noeud1();
441     MG_NOEUD* noeud2=tri->get_noeud2();
442     MG_NOEUD* noeud3=tri->get_noeud3();
443     lsttricontraint.supprimer(tri);
444     delete tri;
445     delete noeud1;
446     delete noeud2;
447     delete noeud3;
448     }
449     }
450     multimap<unsigned long,BLOC_MAILLE_2D*,less <unsigned long> >::iterator itfin2=lstb2d.begin();
451     while (itfin2!=lstb2d.end())
452     {
453     delete (*itfin2).second;
454     itfin2++;
455     }
456     cree_liste_frontiere(maimodifie);
457     // etape 10 Destruction autour des nouveaux triangles
458     if (etape<10) return;
459     affiche(" Destruction autour des nouvelles faces");
460     lstnoeuddetruire.vide();
461     LISTE_MG_TRIANGLE::iterator itnvtri;
462     for (MG_TRIANGLE* tri=maimodifie->get_premier_triangle(itnvtri);tri!=NULL;tri=maimodifie->get_suivant_triangle(itnvtri))
463     {
464     if (tri->get_origine()==MAILLEUR_AUTO)
465     {
466     lstnoeuddetruire.ajouter(tri->get_noeud1());
467     lstnoeuddetruire.ajouter(tri->get_noeud2());
468     lstnoeuddetruire.ajouter(tri->get_noeud3());
469     }
470     }
471     for (MG_NOEUD* no=lstnoeuddetruire.get_premier(it);no!=NULL;no=lstnoeuddetruire.get_suivant(it))
472     {
473     detruit_noeud(no,-1,1);
474     }
475     // etape 11 Creation des blocs de mailles 3D
476     if (etape<11) return;
477     affiche(" Creation des blocs de maille 3D");
478     int nbvolume=geoorigine->get_nb_mg_volume();
479     for (int i=0;i<nbvolume;i++)
480     {
481     multimap<unsigned long,MG_TETRA*,less <unsigned long> > lsttet;
482     MG_VOLUME* vol=geoorigine->get_mg_volume(i);
483     int nblien=vol->get_lien_maillage()->get_nb();
484     TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR ittet;
485     for (MG_ELEMENT_MAILLAGE* ele=vol->get_lien_maillage()->get_premier(ittet);ele!=NULL;ele=vol->get_lien_maillage()->get_suivant(ittet))
486     {
487     MG_TETRA* tet=(MG_TETRA*)ele;
488     if ((maiorigine->get_mg_tetraid(tet->get_id())!=NULL))
489     {
490     MG_NOEUD* n1=tet->get_noeud1();
491     MG_NOEUD* n2=tet->get_noeud2();
492     MG_NOEUD* n3=tet->get_noeud3();
493     MG_NOEUD* n4=tet->get_noeud4();
494     if (n1->get_nouveau_numero()==CONSERVE)
495     if (n2->get_nouveau_numero()==CONSERVE)
496     if (n3->get_nouveau_numero()==CONSERVE)
497     if (n4->get_nouveau_numero()==CONSERVE)
498     {
499     tet->change_nouveau_numero(NONTRAITE);
500     pair<unsigned long,MG_TETRA*> tmp(n1->get_id(),tet);
501     lsttet.insert(tmp);
502     pair<unsigned long,MG_TETRA*> tmp2(n2->get_id(),tet);
503     lsttet.insert(tmp2);
504     pair<unsigned long,MG_TETRA*> tmp3(n3->get_id(),tet);
505     lsttet.insert(tmp3);
506     pair<unsigned long,MG_TETRA*> tmp4(n4->get_id(),tet);
507     lsttet.insert(tmp4);
508     }
509     }
510    
511     }
512     cree_bloc_maille(lsttet,vol);
513     }
514     // etape 12 Maillage 3D
515     if (etape<12) return;
516     affiche(" MAILLEUR 3D");
517     nbvolume=geomodifie->get_nb_mg_volume();
518     for (int i=0;i<nbvolume;i++)
519     {
520     TPL_LISTE_ENTITE<MG_TETRA*> lsttetcontraint;
521     MG_VOLUME* vol=geomodifie->get_mg_volume(i);
522     multimap<unsigned long,BLOC_MAILLE_3D*,less <unsigned long> >::iterator it=lstb3d.begin();
523     while (it!=lstb3d.end())
524     {
525     BLOC_MAILLE_3D *b3d=(*it).second;
526     if (b3d->etat==NONATTACHE)
527     {
528     MG_NOEUD* no1=b3d->lst[0]->get_noeud1();
529     MG_NOEUD* no2=b3d->lst[0]->get_noeud2();
530     MG_NOEUD* no3=b3d->lst[0]->get_noeud3();
531     MG_NOEUD* no4=b3d->lst[0]->get_noeud4();
532     double xyz1[3],xyz2[3],xyz3[3],xyz4[3],xyz[3];
533     transfert_coord(no1->get_coord(),xyz1);
534     transfert_coord(no2->get_coord(),xyz2);
535     transfert_coord(no3->get_coord(),xyz3);
536     transfert_coord(no4->get_coord(),xyz4);
537     xyz[0]=0.25*(xyz1[0]+xyz2[0]+xyz3[0]+xyz4[0]);
538     xyz[1]=0.25*(xyz1[1]+xyz2[1]+xyz3[1]+xyz4[1]);
539     xyz[2]=0.25*(xyz1[2]+xyz2[2]+xyz3[2]+xyz4[2]);
540     int res=point_appartient_volume(xyz,i);
541     if (res==1)
542     {
543     int nbtet=b3d->lst.size();
544     b3d->etat=ATTACHE;
545     for (int j=0;j<nbtet;j++)
546     {
547     MG_TETRA *tet=b3d->lst[j];
548     MG_NOEUD* no1=tet->get_noeud1();
549     MG_NOEUD* no2=tet->get_noeud2();
550     MG_NOEUD* no3=tet->get_noeud3();
551     MG_NOEUD* no4=tet->get_noeud4();
552     double xyz1[3],xy2[3],xyz3[3],xyz4[3];
553     transfert_coord(no1->get_coord(),xyz1);
554     transfert_coord(no2->get_coord(),xyz2);
555     transfert_coord(no3->get_coord(),xyz3);
556     transfert_coord(no4->get_coord(),xyz4);
557     MG_NOEUD* nvno1=new MG_NOEUD(NULL,xyz1[0],xyz1[1],xyz1[2],IMPOSE);
558     MG_NOEUD* nvno2=new MG_NOEUD(NULL,xyz2[0],xyz2[1],xyz2[2],IMPOSE);
559     MG_NOEUD* nvno3=new MG_NOEUD(NULL,xyz3[0],xyz3[1],xyz3[2],IMPOSE);
560     MG_NOEUD* nvno4=new MG_NOEUD(NULL,xyz4[0],xyz4[1],xyz4[2],IMPOSE);
561     MG_TETRA* nvtet=new MG_TETRA(NULL,nvno1,nvno2,nvno3,nvno4,NULL,NULL,NULL,NULL,IMPOSE);
562     lsttetcontraint.ajouter(nvtet);
563     }
564     }
565     }
566     it++;
567     }
568     MAILLEUR3D m3d(maimodifie,geomodifie,vol,cartemod);
569     char mess[50];
570     sprintf(mess," Volume %i identificateur %lu",i,vol->get_id());
571     affiche(mess);
572     m3d.active_affichage(affiche);
573     m3d.maille(vol,NULL,&(lsttetcontraint));
574     int nbtetcon=lsttetcontraint.get_nb();
575     for (int i=0;i<nbtetcon;i++)
576     {
577     MG_TETRA* tet=lsttetcontraint.get(0);
578     MG_NOEUD* noeud1=tet->get_noeud1();
579     MG_NOEUD* noeud2=tet->get_noeud2();
580     MG_NOEUD* noeud3=tet->get_noeud3();
581     MG_NOEUD* noeud4=tet->get_noeud4();
582     lsttetcontraint.supprimer(tet);
583     delete tet;
584     delete noeud1;
585     delete noeud2;
586     delete noeud3;
587     delete noeud4;
588     }
589     }
590     multimap<unsigned long,BLOC_MAILLE_3D*,less <unsigned long> >::iterator itfin3=lstb3d.begin();
591     while (itfin3!=lstb3d.end())
592     {
593     delete (*itfin3).second;
594     itfin3++;
595     }
596    
597    
598     }
599    
600    
601     void REMAILLEUR::initialise_octree(void)
602     {
603     FCT_GENERATEUR_3D<4> *carte2=(FCT_GENERATEUR_3D<4> *)carteori;
604     int nb_cellule=carte2->get_nb_cellule();
605     double param[32];
606     carte2->get_information(0,0,param);
607     double xmin=param[0]+0.005*(param[1]-param[0]);
608     double ymin=param[8]+0.005*(param[10]-param[8]);
609     double zmin=param[16]+0.005*(param[20]-param[16]);
610     carte2->get_information(nb_cellule-1,0,param);
611     double xmax=param[1]-0.005*(param[1]-param[0]);
612     double ymax=param[10]-0.005*(param[10]-param[8]);
613     double zmax=param[20]-0.005*(param[20]-param[16]);
614     octree=new TPL_OCTREE_FCT<MG_NOEUD*,FCT_GENERATEUR_3D<4> >;
615     octree->initialiser(*carte2,xmin,ymin,zmin,xmax,ymax,zmax);
616     LISTE_MG_NOEUD::iterator it;
617     for (MG_NOEUD* noeud=maiorigine->get_premier_noeud(it);noeud!=NULL;noeud=maiorigine->get_suivant_noeud(it))
618     {
619     noeud->change_nouveau_numero(CONSERVE);
620     octree->inserer(noeud);
621     }
622    
623     }
624    
625     void REMAILLEUR::detruit_noeud(MG_NOEUD* noeudbase,double distance,int type)
626     {
627     double* xyz=noeudbase->get_coord();
628     TPL_MAP_ENTITE<MG_NOEUD*> lst;
629     double Ea,Eb;
630     if (type!=0)
631     {
632     double m[9];
633     carteori->evaluer(xyz,m);
634     Ea=2./sqrt(m[0]);
635     distance=2*Ea;
636     }
637     octree->rechercher(xyz[0],xyz[1],xyz[2],distance,lst);
638     LISTE_MG_NOEUD::iterator it;
639     for (MG_NOEUD* noeud=lst.get_premier(it);noeud!=NULL;noeud=lst.get_suivant(it))
640     {
641     double *xyztmp=noeud->get_coord();
642     OT_VECTEUR_3D vec(xyz,xyztmp);
643     if (type!=0)
644     {
645     double m[9];
646     carteori->evaluer(xyztmp,m);
647     Eb=2./sqrt(m[0]);
648     distance=std::max(Ea,Eb);
649     distance=std::max(distance,2*fabs(Ea-Eb));
650     }
651     double cosinus=1.;
652     if (noeudbase->get_lien_topologie()->get_dimension()==2)
653     {
654     MG_FACE* face=(MG_FACE*)noeudbase->get_lien_topologie();
655     double uv[2];
656     double n[3];
657     face->inverser(uv,xyz);
658     face->calcul_normale_unitaire(uv,n);
659     OT_VECTEUR_3D norm(n);
660     cosinus=norm*vec;
661     cosinus=cosinus/vec.get_longueur();
662     cosinus=(1.-0.3)*cosinus+0.3; //attention c est 0.5 dans la these mais 0.3 si ca marche c est mieux
663     }
664     if (vec.get_longueur()<distance*cosinus)
665     noeud->change_nouveau_numero(DETRUIT);
666     }
667    
668     }
669    
670    
671     void REMAILLEUR::cree_quadtree(MG_MAILLAGE* mai,TPL_QUADTREE<MG_SEGMENT_FRONTIERE*,MG_NOEUD*> **quad,TPL_MAP_ENTITE<MG_SEGMENT_FRONTIERE*> *lstsegfront)
672     {
673     int nbface=mai->get_mg_geometrie()->get_nb_mg_face();
674     for (int i=0;i<nbface;i++)
675     {
676     ((quad[i]))=new TPL_QUADTREE<MG_SEGMENT_FRONTIERE*,MG_NOEUD*>;
677     TPL_LISTE_ENTITE<MG_NOEUD*> lstn;
678     double umax=-1e300;
679     double vmax=-1e300;
680     double umin=1e300;
681     double vmin=1e300;
682     MG_FACE* face=mai->get_mg_geometrie()->get_mg_face(i);
683     int nbboucle=face->get_nb_mg_boucle();
684     int nbsegfravant=lstsegfront->get_nb();
685     for (int j=0;j<nbboucle;j++)
686     {
687     MG_BOUCLE* bou=face->get_mg_boucle(j);
688     int nbarete=bou->get_nb_mg_coarete();
689     for (int k=0;k<nbarete;k++)
690     {
691     MG_ARETE* are=bou->get_mg_coarete(k)->get_arete();
692     int nbseg=are->get_lien_maillage()->get_nb();
693     for (int l=0;l<nbseg;l++)
694     {
695     MG_SEGMENT* seg=(MG_SEGMENT*)are->get_lien_maillage()->get(l);
696     if (mai->get_mg_segmentid(seg->get_id())==NULL) continue;
697     double *xyz=seg->get_noeud1()->get_coord();
698     double uv[2];
699     face->inverser(uv,xyz);
700     seg->get_noeud1()->change_u(uv[0]);
701     seg->get_noeud1()->change_v(uv[1]);
702     xyz=seg->get_noeud2()->get_coord();
703     face->inverser(uv,xyz);
704     seg->get_noeud2()->change_u(uv[0]);
705     seg->get_noeud2()->change_v(uv[1]);
706     lstn.ajouter(seg->get_noeud2());
707     if (uv[0]<umin) umin=uv[0];
708     if (uv[0]>umax) umax=uv[0];
709     if (uv[1]<vmin) vmin=uv[1];
710     if (uv[1]>vmax) vmax=uv[1];
711     MG_SEGMENT_FRONTIERE* segfr=new MG_SEGMENT_FRONTIERE(seg);
712     lstsegfront->ajouter(segfr);
713     }
714     }
715     }
716     double periodeenu=face->get_surface()->get_periode_u();
717     double periodeenv=face->get_surface()->get_periode_v();
718     if (periodeenu!=0.0)
719     {
720     umin=0.;
721     umax=periodeenu;
722     }
723     else
724     {
725     double diff=umax-umin;
726     umin=umin-0.125*diff;
727     umax=umax+0.125*diff;
728     }
729     if (periodeenv!=0.0)
730     {
731     vmin=0.;
732     vmax=periodeenv;
733     }
734     else
735     {
736     double diff=vmax-vmin;
737     vmin=vmin-0.125*diff;
738     vmax=vmax+0.125*diff;
739     }
740     (quad[i])->initialiser(&lstn,1,umin,vmin,umax,vmax,face->get_surface()->get_periode_u(),face->get_surface()->get_periode_v());
741     int nbsegfr=lstsegfront->get_nb();
742     for (int j=nbsegfravant;j<nbsegfr;j++)
743     {
744     double *uv=lstsegfront->get(j)->get_uv1();
745     uv[0]=lstsegfront->get(j)->get_segment()->get_noeud1()->get_u();
746     uv[1]=lstsegfront->get(j)->get_segment()->get_noeud1()->get_v();
747     uv=lstsegfront->get(j)->get_uv2();
748     uv[0]=lstsegfront->get(j)->get_segment()->get_noeud2()->get_u();
749     uv[1]=lstsegfront->get(j)->get_segment()->get_noeud2()->get_v();
750     (quad[i])->inserer( lstsegfront->get(j));
751     }
752     }
753    
754     }
755    
756     void REMAILLEUR::cree_liste_frontiere(MG_MAILLAGE* mai)
757     {
758     int nbvolume=mai->get_mg_geometrie()->get_nb_mg_volume();
759     for (int i=0;i<nbvolume;i++)
760     {
761     MG_VOLUME* vol=mai->get_mg_geometrie()->get_mg_volume(i);
762     int nbco=vol->get_nb_mg_coquille();
763     for (int j=0;j<nbco;j++)
764     {
765     MG_COQUILLE* coq=vol->get_mg_coquille(j);
766     int nbface=coq->get_nb_mg_coface();
767     for (int k=0;k<nbface;k++)
768     {
769     MG_FACE* face=coq->get_mg_coface(k)->get_face();
770     int nbseg=face->get_lien_maillage()->get_nb();
771     for (int l=0;l<nbseg;l++)
772     {
773     MG_TRIANGLE* tri=(MG_TRIANGLE*)face->get_lien_maillage()->get(l);
774     if (mai->get_mg_triangleid(tri->get_id())==NULL) continue;
775     lsttrifront[i].ajouter(tri);
776     }
777     }
778     }
779    
780    
781     }
782     }
783    
784    
785    
786     double REMAILLEUR::calcule_distance_contour_face(double *uv,MG_FACE* face,int num,TPL_QUADTREE<MG_SEGMENT_FRONTIERE*,MG_NOEUD*> **quad)
787     {
788     double uvajoute[3]={uv[0],uv[1],0.};
789     TPL_MAP_ENTITE<MG_COARETE*> lst;
790     int nbboucle=face->get_nb_mg_boucle();
791     for (int i=0;i<nbboucle;i++)
792     {
793     int nbcoarete=face->get_mg_boucle(i)->get_nb_mg_coarete();
794     for (int j=0;j<nbcoarete;j++)
795     lst.ajouter(face->get_mg_boucle(i)->get_mg_coarete(j));
796     }
797     BOITE_2D b=quad[num]->get_cellule(0)->get_boite();
798     double rayon=b.get_rayon()/quad[num]->get_nb_cellule();
799     TPL_MAP_ENTITE<MG_SEGMENT_FRONTIERE*> lsttrouve;
800     MG_SEGMENT_FRONTIERE* segretenu=NULL;
801     MG_COARETE* coareteretenu=NULL;
802     double distance=1e300;
803     do
804     {
805     do
806     {
807     quad[num]->rechercher(uv[0],uv[1],rayon,lsttrouve);
808     rayon=rayon*1.25;
809     }
810     while (lsttrouve.get_nb()==0);
811     int nb=lsttrouve.get_nb();
812     for (int i=0;i<nb;i++)
813     {
814     MG_SEGMENT_FRONTIERE* seg=lsttrouve.get(i);
815     MG_ARETE* are=(MG_ARETE*)seg->get_segment()->get_lien_topologie();
816     int nbco=are->get_nb_mg_coarete();
817     MG_COARETE* coarete=NULL;
818     for (int j=0;j<nbco;j++)
819     if (lst.getid(are->get_mg_coarete(j)->get_id())!=NULL) coarete=are->get_mg_coarete(j);
820     if (coarete!=NULL)
821     {
822     double uv1[3],uv2[3];
823     uv1[0]=seg->get_uv1()[0];
824     uv1[1]=seg->get_uv1()[1];
825     uv1[2]=0;
826     uv2[0]=seg->get_uv2()[0];
827     uv2[1]=seg->get_uv2()[1];
828     uv2[2]=0;
829     double dist=distance_pt_segment(uvajoute,uv1,uv2,face);
830     if (dist<distance)
831     {
832     distance=dist;
833     segretenu=seg;
834     coareteretenu=coarete;
835     }
836     }
837     }
838     }
839     while (rayon<distance);
840     double uv1[3]={segretenu->get_uv1()[0],segretenu->get_uv1()[1],0.};
841     double uv2[3]={segretenu->get_uv2()[0],segretenu->get_uv2()[1],0.};
842     OT_DECALAGE_PARAMETRE deca(face->get_surface()->get_periode_u(),face->get_surface()->get_periode_v());
843     double du=deca.calcul_decalage_parametre_u(uv1[0]);
844     double dv=deca.calcul_decalage_parametre_v(uv1[1]);
845     uv1[0]=uv1[0]+du;
846     uv1[1]=uv1[1]+dv;
847     uv2[0]=deca.decalage_parametre_u(uv2[0],du);
848     uv2[1]=deca.decalage_parametre_v(uv2[1],dv);
849     uvajoute[0]=deca.decalage_parametre_u(uvajoute[0],du);
850     uvajoute[1]=deca.decalage_parametre_v(uvajoute[1],dv);
851     OT_VECTEUR_3D base(uv1,uv2);
852     base=coareteretenu->get_orientation()*base;
853     double normal[3]={0.,0.,1.};
854     OT_VECTEUR_3D nor(normal);
855     OT_VECTEUR_3D dir(uv1,uvajoute);
856     base.norme();
857     dir.norme();
858     double sens=(nor&base)*dir;
859     if (sens<0.) distance=-distance;
860     return distance;
861     }
862    
863    
864     double REMAILLEUR::distance_pt_segment(double *uv,MG_NOEUD* noeud1,MG_NOEUD *noeud2,MG_FACE* face)
865     {
866     double uv1[3],uv2[3];
867     face->inverser(uv1,noeud1->get_coord());
868     noeud1->change_u(uv1[0]);
869     noeud1->change_v(uv1[1]);
870     face->inverser(uv2,noeud2->get_coord());
871     noeud2->change_u(uv2[0]);
872     noeud2->change_v(uv2[1]);
873     uv1[2]=0.;uv2[2]=0.;
874     return distance_pt_segment(uv,uv1,uv2,face);
875     }
876    
877     double REMAILLEUR::distance_pt_segment(double *uv,double *uv1,double *uv2,MG_FACE* face)
878     {
879     OT_DECALAGE_PARAMETRE deca(face->get_surface()->get_periode_u(),face->get_surface()->get_periode_v());
880     double du=deca.calcul_decalage_parametre_u(uv1[0]);
881     double dv=deca.calcul_decalage_parametre_v(uv1[1]);
882     uv1[0]=uv1[0]+du;
883     uv1[1]=uv1[1]+dv;
884     uv2[0]=deca.decalage_parametre_u(uv2[0],du);
885     uv2[1]=deca.decalage_parametre_v(uv2[1],dv);
886     double uvtmp[3];
887     uvtmp[0]=deca.decalage_parametre_u(uv[0],du);
888     uvtmp[1]=deca.decalage_parametre_v(uv[1],dv);
889     uvtmp[2]=0.;
890     double distance;
891     OT_VECTEUR_3D vec12(uv1,uv2);
892     OT_VECTEUR_3D vec1p(uv1,uvtmp);
893     double t=(vec12.get_x()*vec1p.get_x()+vec12.get_y()*vec1p.get_y()+vec12.get_z()*vec1p.get_z())/(vec12.get_x()*vec12.get_x()+vec12.get_y()*vec12.get_y()+vec12.get_z()*vec12.get_z());
894     if ((t<0.) || (t>1.))
895     {
896     OT_VECTEUR_3D vec2p(uv2,uvtmp);
897     double dist1=vec1p.get_longueur();
898     double dist2=vec2p.get_longueur();
899     distance=0.5*(dist1+dist2);
900     }
901     else
902     {
903     OT_VECTEUR_3D pvec=vec12&vec1p;
904     distance=pvec.get_longueur()/vec12.get_longueur();
905     }
906     return distance;
907     }
908    
909     void REMAILLEUR::ajouter_liste_correspondance_noeud(MG_NOEUD* n1,MG_NOEUD* n2)
910     {
911     CORRESPONDANCENOEUD corr(n1,n2);
912     pair<unsigned long,CORRESPONDANCENOEUD> tmp(n1->get_id(),corr);
913     lstcorrnoeud.insert(tmp);
914     }
915    
916     MG_NOEUD* REMAILLEUR::get_liste_correspondance(MG_NOEUD* n1)
917     {
918     MG_NOEUD *n2=NULL;
919     std::map<unsigned long,CORRESPONDANCENOEUD,std::less<unsigned long> >::iterator j=lstcorrnoeud.find(n1->get_id());
920     if (j!=lstcorrnoeud.end())
921     n2=(*j).second.noeudmod;
922     return n2;
923     }
924    
925    
926     MG_NOEUD* REMAILLEUR::transfert_noeud(MG_NOEUD* no,MG_ELEMENT_TOPOLOGIQUE* ele)
927     {
928     double *xyz=no->get_coord();
929     OT_VECTEUR_3D pt(xyz);
930     OT_VECTEUR_3D nvpt=rescmp.get_coord(pt);
931     MG_NOEUD* noeud=new MG_NOEUD(ele,nvpt.get_x(),nvpt.get_y(),nvpt.get_z(),IMPOSE);
932     maimodifie->ajouter_mg_noeud(noeud);
933     ajouter_liste_correspondance_noeud(no,noeud);
934     return noeud;
935     }
936    
937     void REMAILLEUR::transfert_coord(double *xyz,double *xyz2)
938     {
939     OT_VECTEUR_3D pt(xyz);
940     OT_VECTEUR_3D nvpt=rescmp.get_coord(pt);
941     xyz2[0]=nvpt.get_x();
942     xyz2[1]=nvpt.get_y();
943     xyz2[2]=nvpt.get_z();
944     return;
945     }
946    
947     void REMAILLEUR::cree_bloc_maille(multimap<unsigned long,MG_SEGMENT*,less <unsigned long> > &lst,MG_ARETE* arete)
948     {
949     while (lst.size()>0)
950     {
951     multimap<unsigned long,MG_SEGMENT*,less <unsigned long> >::iterator it=lst.begin();
952     MG_SEGMENT* seg=(*it).second;
953     lst.erase(it);
954     TPL_LISTE_ENTITE<double> param;
955     BLOC_MAILLE_1D *b1d=new BLOC_MAILLE_1D(arete,arete->get_courbe()->get_type_geometrique(param));
956     pair<unsigned long,BLOC_MAILLE_1D*> tmp(b1d->get_id(),b1d);
957     lstb1d.insert(tmp);
958     b1d->lst.push_back(seg);
959     seg->change_nouveau_numero(TRAITE);
960     //char mess[255];
961     //sprintf(mess," Arete %lu , Bloc %d",arete->get_id(),lstb1d.size());
962     //affiche(mess);
963     int num=0;
964     while (num!=b1d->lst.size())
965     {
966     MG_SEGMENT* itseg=b1d->lst[num++];
967     MG_SEGMENT* seg=trouve_segment(lst,itseg->get_noeud1());
968     while (seg!=NULL)
969     {
970     if (seg->get_nouveau_numero()==NONTRAITE)
971     {
972     b1d->lst.push_back(seg);
973     seg->change_nouveau_numero(TRAITE);
974     }
975     seg=trouve_segment(lst,itseg->get_noeud1());
976     }
977     seg=trouve_segment(lst,itseg->get_noeud2());
978     while (seg!=NULL)
979     {
980     if (seg->get_nouveau_numero()==NONTRAITE)
981     {
982     b1d->lst.push_back(seg);
983     seg->change_nouveau_numero(TRAITE);
984     }
985     seg=trouve_segment(lst,itseg->get_noeud2());
986     }
987     }
988    
989     //sprintf(mess," Segment %d ",b1d->lst.size());
990     //affiche(mess);
991     }
992     }
993    
994     void REMAILLEUR::cree_bloc_maille(multimap<unsigned long,MG_TRIANGLE*,less <unsigned long> > &lst,MG_FACE* face)
995     {
996     while (lst.size()>0)
997     {
998     multimap<unsigned long,MG_TRIANGLE*,less <unsigned long> >::iterator it=lst.begin();
999     MG_TRIANGLE* tri=(*it).second;
1000     lst.erase(it);
1001     TPL_LISTE_ENTITE<double> param;
1002     BLOC_MAILLE_2D *b2d=new BLOC_MAILLE_2D(geoorigine->get_mg_faceid(face->get_id()),face->get_surface()->get_type_geometrique(param));
1003     pair<unsigned long,BLOC_MAILLE_2D*> tmp(b2d->get_id(),b2d);
1004     lstb2d.insert(tmp);
1005     b2d->lst.push_back(tri);
1006     tri->change_nouveau_numero(TRAITE);
1007     //char mess[255];
1008     //sprintf(mess," Face %lu , Bloc %d",face->get_id(),lstb2d.size());
1009     //affiche(mess);
1010     int num=0;
1011     while (num!=b2d->lst.size())
1012     {
1013     MG_TRIANGLE* ittri=b2d->lst[num++];
1014     MG_TRIANGLE* tri=trouve_triangle(lst,ittri->get_noeud1());
1015     while (tri!=NULL)
1016     {
1017     if (tri->get_nouveau_numero()==NONTRAITE)
1018     {
1019     b2d->lst.push_back(tri);
1020     tri->change_nouveau_numero(TRAITE);
1021     }
1022     tri=trouve_triangle(lst,ittri->get_noeud1());
1023     }
1024     tri=trouve_triangle(lst,ittri->get_noeud2());
1025     while (tri!=NULL)
1026     {
1027     if (tri->get_nouveau_numero()==NONTRAITE)
1028     {
1029     b2d->lst.push_back(tri);
1030     tri->change_nouveau_numero(TRAITE);
1031     }
1032     tri=trouve_triangle(lst,ittri->get_noeud2());
1033     }
1034     tri=trouve_triangle(lst,ittri->get_noeud3());
1035     while (tri!=NULL)
1036     {
1037     if (tri->get_nouveau_numero()==NONTRAITE)
1038     {
1039     b2d->lst.push_back(tri);
1040     tri->change_nouveau_numero(TRAITE);
1041     }
1042     tri=trouve_triangle(lst,ittri->get_noeud3());
1043     }
1044     }
1045    
1046     //sprintf(mess," Triangle %d ",b2d->lst.size());
1047     //affiche(mess);
1048     }
1049     }
1050    
1051     void REMAILLEUR::cree_bloc_maille(multimap<unsigned long,MG_TETRA*,less <unsigned long> > &lst,MG_VOLUME* vol)
1052     {
1053     while (lst.size()>0)
1054     {
1055     multimap<unsigned long,MG_TETRA*,less <unsigned long> >::iterator it=lst.begin();
1056     MG_TETRA* tet=(*it).second;
1057     lst.erase(it);
1058     BLOC_MAILLE_3D *b3d=new BLOC_MAILLE_3D(geoorigine->get_mg_volumeid(vol->get_id()));
1059     pair<unsigned long,BLOC_MAILLE_3D*> tmp(b3d->get_id(),b3d);
1060     lstb3d.insert(tmp);
1061     b3d->lst.push_back(tet);
1062     tet->change_nouveau_numero(TRAITE);
1063     //char mess[255];
1064     //sprintf(mess," Volume %lu , Bloc %d",vol->get_id(),lstb3d.size());
1065     //affiche(mess);
1066     int num=0;
1067     while (num!=b3d->lst.size())
1068     {
1069     MG_TETRA* ittet=b3d->lst[num++];
1070     MG_TETRA* tet=trouve_tetra(lst,ittet->get_noeud1());
1071     while (tet!=NULL)
1072     {
1073     if (tet->get_nouveau_numero()==NONTRAITE)
1074     {
1075     b3d->lst.push_back(tet);
1076     tet->change_nouveau_numero(TRAITE);
1077     }
1078     tet=trouve_tetra(lst,ittet->get_noeud1());
1079     }
1080     tet=trouve_tetra(lst,ittet->get_noeud2());
1081     while (tet!=NULL)
1082     {
1083     if (tet->get_nouveau_numero()==NONTRAITE)
1084     {
1085     b3d->lst.push_back(tet);
1086     tet->change_nouveau_numero(TRAITE);
1087     }
1088     tet=trouve_tetra(lst,ittet->get_noeud2());
1089     }
1090     tet=trouve_tetra(lst,ittet->get_noeud3());
1091     while (tet!=NULL)
1092     {
1093     if (tet->get_nouveau_numero()==NONTRAITE)
1094     {
1095     b3d->lst.push_back(tet);
1096     tet->change_nouveau_numero(TRAITE);
1097     }
1098     tet=trouve_tetra(lst,ittet->get_noeud3());
1099     }
1100     tet=trouve_tetra(lst,ittet->get_noeud4());
1101     while (tet!=NULL)
1102     {
1103     if (tet->get_nouveau_numero()==NONTRAITE)
1104     {
1105     b3d->lst.push_back(tet);
1106     tet->change_nouveau_numero(TRAITE);
1107     }
1108     tet=trouve_tetra(lst,ittet->get_noeud4());
1109     }
1110     }
1111    
1112     //sprintf(mess," Tetra %d ",b3d->lst.size());
1113     //affiche(mess);
1114     }
1115     }
1116    
1117    
1118     MG_SEGMENT* REMAILLEUR::trouve_segment(multimap<unsigned long,MG_SEGMENT*,less <unsigned long> > &lst,MG_NOEUD* no)
1119     {
1120     MG_SEGMENT* seg=NULL;
1121     unsigned long id=no->get_id();
1122     multimap<unsigned long,MG_SEGMENT*,less <unsigned long> >::iterator it=lst.find(id);
1123     if (it!=lst.end())
1124     {
1125     seg=(*it).second;
1126     lst.erase(it);
1127     }
1128     return seg;
1129     }
1130    
1131     MG_TRIANGLE* REMAILLEUR::trouve_triangle(multimap<unsigned long,MG_TRIANGLE*,less <unsigned long> > &lst,MG_NOEUD* no)
1132     {
1133     MG_TRIANGLE* tri=NULL;
1134     unsigned long id=no->get_id();
1135     multimap<unsigned long,MG_TRIANGLE*,less <unsigned long> >::iterator it=lst.find(id);
1136     if (it!=lst.end())
1137     {
1138     tri=(*it).second;
1139     lst.erase(it);
1140     }
1141     return tri;
1142     }
1143    
1144     MG_TETRA* REMAILLEUR::trouve_tetra(multimap<unsigned long,MG_TETRA*,less <unsigned long> > &lst,MG_NOEUD* no)
1145     {
1146     MG_TETRA* tet=NULL;
1147     unsigned long id=no->get_id();
1148     multimap<unsigned long,MG_TETRA*,less <unsigned long> >::iterator it=lst.find(id);
1149     if (it!=lst.end())
1150     {
1151     tet=(*it).second;
1152     lst.erase(it);
1153     }
1154     return tet;
1155     }
1156    
1157     int REMAILLEUR::point_appartient_volume(double *xyz,int numvol)
1158     {
1159     int ok=0;
1160     int i=0;
1161     multimap<double,double,less<double> > intersection;
1162     OT_VECTEUR_3D vecteur_dir;
1163     while (ok==0)
1164    
1165     {
1166     int ok2=0;
1167     while (ok2==0)
1168     {
1169     MG_TRIANGLE* tri=lsttrifront[numvol].get(i);
1170     MG_NOEUD* no1=tri->get_noeud1();
1171     MG_NOEUD* no2=tri->get_noeud2();
1172     MG_NOEUD* no3=tri->get_noeud3();
1173     double xb=(no1->get_x()+no2->get_x()+no3->get_x())/3.;
1174     double yb=(no1->get_y()+no2->get_y()+no3->get_y())/3.;
1175     double zb=(no1->get_z()+no2->get_z()+no3->get_z())/3.;
1176     OT_VECTEUR_3D directeur(xb-xyz[0],yb-xyz[1],zb-xyz[2]);
1177     directeur.norme();
1178     vecteur_dir=directeur;
1179     OT_VECTEUR_3D n1n3(no1->get_coord(),no3->get_coord());
1180     OT_VECTEUR_3D n1n2(no1->get_coord(),no2->get_coord());
1181     OT_VECTEUR_3D n=n1n3&n1n2;
1182     n.norme();
1183     double ps1=n*directeur;
1184     if (fabs(ps1)<1e-6) i++;
1185     else ok2=1;
1186     }
1187     ok2=0;
1188     intersection.clear();
1189     for (int j=0;j<lsttrifront[numvol].get_nb();j++)
1190     {
1191     MG_TRIANGLE* tri=lsttrifront[numvol].get(j);
1192    
1193     double t;
1194     int res=inter_droite_triangle(xyz,vecteur_dir,tri,&t);
1195     if (res==2) {ok2=1;break;}
1196     if (res==1)
1197     {
1198     pair<double,double> tmp(t,t);
1199     intersection.insert(tmp);
1200     }
1201     }
1202     if (ok2==0) ok=1;
1203     i++;
1204    
1205     }
1206     multimap<double,double,less <double> >::iterator it=intersection.begin();
1207     int nbinterneg=0;
1208     while (it!=intersection.end())
1209     {
1210     if ((*it).second<0.) nbinterneg++;
1211     it++;
1212     }
1213     if (nbinterneg%2==1) return 1;
1214     return 0;
1215     }
1216    
1217    
1218     int REMAILLEUR::inter_droite_triangle(double *xyz,double *dir,MG_TRIANGLE* tri,double *t)
1219     {
1220     int inter=0;
1221     MG_NOEUD *noeud1=tri->get_noeud1();
1222     MG_NOEUD *noeud2=tri->get_noeud2();
1223     MG_NOEUD *noeud3=tri->get_noeud3();
1224     OT_VECTEUR_3D g1(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
1225     OT_VECTEUR_3D g2(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
1226     OT_VECTEUR_3D g3(dir);
1227     OT_MATRICE_3D systeme(g1,g2,g3);
1228     double det=systeme.get_determinant();
1229     double eps=0.333333333333*(g1.diff()+g2.diff()+g3.diff());
1230     eps=eps*eps*0.0018;
1231     if (OPERATEUR::egal(det,0.0,eps)==true)
1232     {
1233     OT_VECTEUR_3D g3b(xyz[0]-noeud1->get_x(),xyz[1]-noeud1->get_y(),xyz[2]-noeud1->get_z());
1234     OT_MATRICE_3D systeme2(g1,g2,g3b);
1235     double det2=systeme2.get_determinant();
1236     double eps2=0.333333333333*(g1.diff()+g2.diff()+g3b.diff());
1237     eps2=18.*eps2*eps2*eps2*1e-6;
1238     if (OPERATEUR::egal(det2,0.0,eps)==true) // cas 2D
1239     {
1240     return 2;
1241     }
1242     else return 0;
1243     }
1244     else
1245     {
1246     double x,y,z;
1247     x=1.0/det*(g2.get_y()*g3.get_z()-g2.get_z()*g3.get_y());
1248     y=1.0/det*(g3.get_x()*g2.get_z()-g2.get_x()*g3.get_z());
1249     z=1.0/det*(g2.get_x()*g3.get_y()-g2.get_y()*g3.get_x());
1250     OT_VECTEUR_3D g1b(x,y,z);
1251     x=1.0/det*(g3.get_y()*g1.get_z()-g1.get_y()*g3.get_z());
1252     y=1.0/det*(g1.get_x()*g3.get_z()-g3.get_x()*g1.get_z());
1253     z=1.0/det*(g3.get_x()*g1.get_y()-g1.get_x()*g3.get_y());
1254     OT_VECTEUR_3D g2b(x,y,z);
1255     x=1.0/det*(g1.get_y()*g2.get_z()-g1.get_z()*g2.get_y());
1256     y=1.0/det*(g2.get_x()*g1.get_z()-g1.get_x()*g2.get_z());
1257     z=1.0/det*(g1.get_x()*g2.get_y()-g1.get_y()*g2.get_x());
1258     OT_VECTEUR_3D g3b(x,y,z);
1259     OT_VECTEUR_3D n1n4(xyz[0]-noeud1->get_x(),xyz[1]-noeud1->get_y(),xyz[2]-noeud1->get_z());
1260     double alpha1=n1n4*g1b;
1261     double alpha2=n1n4*g2b;
1262     double alpha3=-(n1n4*g3b);
1263     double alpha4=1-alpha1-alpha2;
1264     double eps=0.000001;
1265     if ((alpha1>-eps) && (alpha1<1.+eps))
1266     if ((alpha2>-eps) && (alpha2<1.+eps))
1267     if ((alpha4>-eps) && (alpha4<1.+eps)) {*t=alpha3;return 1;}
1268     return 0;
1269     }
1270    
1271     }