ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/remailleur.cpp
Revision: 155
Committed: Fri Sep 12 14:48:21 2008 UTC (16 years, 8 months ago) by francois
Original Path: magic/lib/mailleur/mailleur/src/remailleur.cpp
File size: 46802 byte(s)
Log Message:
Mise a jour remailleur il y a encore un pb avec la distance de destruction qui fait que le remaillage ne donne pas des qualite top

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