ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_bloc.cpp
Revision: 120
Committed: Tue Jul 8 19:31:35 2008 UTC (16 years, 10 months ago) by francois
Original Path: magic/lib/mailleur/mailleur/src/mailleur_bloc.cpp
File size: 24795 byte(s)
Log Message:
ajout du mailleur par bloc pour l optimisation

File Contents

# User Rev Content
1 francois 120 //---------------------------------------------------------------------------
2    
3     #include "gestionversion.h"
4     #include <math.h>
5     #include "mailleur_bloc.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     MAILLEUR_BLOC::MAILLEUR_BLOC(class MG_GESTIONNAIRE* gentier,MG_GESTIONNAIRE* gbloc,int numgeoentier,int numgeobloc,class FCT_TAILLE* fct_taille):gestentier(gentier),gestbloc(gbloc),carte(fct_taille)
20     {
21     geoentier=gestentier->get_mg_geometrie(numgeoentier);
22     geobloc=gestbloc->get_mg_geometrie(numgeobloc);
23     }
24    
25     MAILLEUR_BLOC::~MAILLEUR_BLOC()
26     {
27     delete octreeentierseg;
28     delete octreeblocseg;
29     delete octreeentier;
30     delete octreebloc;
31     }
32    
33    
34     void MAILLEUR_BLOC::maille(void)
35     {
36     //maillage 1d temporaire pour fonctionnement et initialisation
37     maitemp=new MG_MAILLAGE(geobloc);
38     maibloc=new MG_MAILLAGE(geobloc);
39     maientier=new MG_MAILLAGE(geoentier);
40     gestbloc->ajouter_mg_maillage(maitemp);
41     gestbloc->ajouter_mg_maillage(maibloc);
42     gestentier->ajouter_mg_maillage(maientier);
43     MAILLEUR0D m0d(maitemp,geobloc);
44     m0d.maille();
45     MAILLEUR1D m1d(maitemp,geobloc,NULL,carte);
46     m1d.maille();
47     cree_octree(maitemp,&octreeblocseg,&octreebloc);
48     // etape 1 maillage 0d du modele entier
49     affiche(" MAILLAGE 0D du modèle entier");
50     MAILLEUR0D m0d_entier(maientier,geoentier);
51     m0d_entier.maille();
52     // etape 2 transfert noeud vers bloc
53     //int nbaretebloc=geobloc->get_nb_mg_arete();
54     //LISTE_NOEUD_ADD_ARETE *lst_noeuds_bloc=new LISTE_NOEUD_ADD_ARETE[nbaretebloc];
55     int nbnoeudentier=maientier->get_nb_mg_noeud();
56     for (int i=0;i<nbnoeudentier;i++)
57     {
58     MG_NOEUD* noeud=maientier->get_mg_noeud(i);
59     double *xyz=noeud->get_coord();
60     MG_SOMMET* som=verif_noeud_sur_sommet(xyz, geobloc);
61     if (som!=NULL)
62     {
63     MG_NOEUD *nvnoeud=new MG_NOEUD(som, xyz[0],xyz[1],xyz[2], IMPOSE);
64     maibloc->ajouter_mg_noeud(nvnoeud);
65     octreebloc->inserer(nvnoeud);
66     }
67     /* else
68     {
69     int num;
70     double t;
71     MG_ARETE* are=verif_noeud_sur_arete(xyz, geobloc, &num,&t);
72     if (are!=NULL)
73     {
74     MG_NOEUD_ADDITIONNEL_ARETE noeud_add;
75     noeud_add.noeud=noeud;
76     noeud_add.t=t;
77     std::pair<double,MG_NOEUD_ADDITIONNEL_ARETE> tmp(t,noeud_add);
78     lst_noeuds_bloc[num].insert(tmp);
79     }
80     }*/
81    
82    
83     }
84     // etape 3 maillage 0d du modele bloc
85     affiche(" MAILLAGE 0D du modele bloc");
86     int nbnoeudblocavantm0d=maibloc->get_nb_mg_noeud();
87     MAILLEUR0D m0d_bloc(maibloc,geobloc);
88     m0d_bloc.adapte();
89     // etape 4 transfert noeud vers entier
90     int nbnouveaunoeudbloc=maibloc->get_nb_mg_noeud()-nbnoeudblocavantm0d;
91     int nbareteentier=geoentier->get_nb_mg_arete();
92     LISTE_NOEUD_ADD_ARETE *lst_noeuds_entier=new LISTE_NOEUD_ADD_ARETE[nbareteentier];
93     for (int i=0;i<nbnouveaunoeudbloc;i++)
94     {
95     MG_NOEUD *noeud=maibloc->get_mg_noeud(i+nbnoeudblocavantm0d);
96     octreebloc->inserer(noeud);
97     double *xyz=noeud->get_coord();
98     int num;double t;
99     MG_ARETE* are=verif_noeud_sur_arete(xyz, geoentier, &num,&t);
100     if (are!=NULL)
101     {
102     MG_NOEUD_ADDITIONNEL_ARETE noeud_add;
103     noeud_add.noeud=noeud;
104     noeud_add.t=t;
105     pair<double,MG_NOEUD_ADDITIONNEL_ARETE > tmp(t,noeud_add);
106     lst_noeuds_entier[num].insert(tmp);
107     }
108     }
109     // etape 5 Maillage 1D de entier
110     affiche(" MAILLAGE 1D du modèle entier");
111     for (int i=0;i<nbareteentier;i++)
112     {
113     MG_ARETE* are=geoentier->get_mg_arete(i);
114     int nbnoeudcontraint=lst_noeuds_entier[i].size();
115     MG_NOEUD* noeuddepart;
116     int numsom=are->get_cosommet1()->get_sommet()->get_lien_maillage()->get_nb();
117     for (int j=0;j<numsom;j++)
118     if (maientier->get_mg_noeudid(are->get_cosommet1()->get_sommet()->get_lien_maillage()->get(j)->get_id())!=NULL) noeuddepart=(MG_NOEUD*)are->get_cosommet1()->get_sommet()->get_lien_maillage()->get(j);
119     MG_NOEUD* noeudarrivee;
120     numsom=are->get_cosommet2()->get_sommet()->get_lien_maillage()->get_nb();
121     for (int j=0;j<numsom;j++)
122     if (maientier->get_mg_noeudid(are->get_cosommet2()->get_sommet()->get_lien_maillage()->get(j)->get_id())!=NULL) noeudarrivee=(MG_NOEUD*)are->get_cosommet2()->get_sommet()->get_lien_maillage()->get(j);
123     MG_NOEUD* navant;
124     double tavant;
125     IT_LISTE_NOEUD_ADD_ARETE it=lst_noeuds_entier[i].begin();
126     for (int j=0;j<nbnoeudcontraint+1;j++)
127     {
128     MG_NOEUD* n1,*n2;
129     double t1,t2;
130     if (j==0) {n1=noeuddepart;t1=are->get_tmin();} else {n1=navant;t1=tavant;}
131     if (j==nbnoeudcontraint) {n2=noeudarrivee;t2=are->get_tmax();}
132     else
133     {
134     double *xyz=(*it).second.noeud->get_coord();
135     n2=new MG_NOEUD(are, xyz[0],xyz[1],xyz[2], IMPOSE);
136     maientier->ajouter_mg_noeud(n2);
137     t2=(*it).second.t;
138     it++;
139     }
140     MAILLEUR1D m1d(maientier, geoentier, are, carte);
141     m1d.maille(are,t1,n1,t2,n2);
142     tavant=t2;
143     navant=n2;
144     }
145     }
146     cree_octree(maientier,&octreeentierseg,&octreeentier);
147     int nbnoeudtmp=maientier->get_nb_mg_noeud();
148     for (int i=0;i<nbnoeudtmp;i++)
149     octreeentier->inserer(maientier->get_mg_noeud(i));
150     //etape 6 transfert des segments vers le maillage bloc
151     //int nbfacebloc=geobloc->get_nb_mg_face();
152     //LISTE_SEGMENT_ADD_FACE *lstsegmentbloc=new LISTE_SEGMENT_ADD_FACE[nbfacebloc];
153     int nbsegment=maientier->get_nb_mg_segment();
154     for (int i=0;i<nbsegment;i++)
155     {
156     MG_SEGMENT* seg=maientier->get_mg_segment(i);
157     MG_NOEUD* noeud1=seg->get_noeud1();
158     MG_NOEUD* noeud2=seg->get_noeud2();
159     int dim1;
160     int dim2;
161     MG_NOEUD* n1=chercher(noeud1->get_coord(),*octreebloc);
162     MG_NOEUD* n2=chercher(noeud2->get_coord(),*octreebloc);
163     MG_ARETE* are1ref;
164     MG_FACE* face1ref;
165     MG_ARETE* are2ref;
166     MG_FACE* face2ref;
167     int num1,num2;
168     double t1,t2;
169     if (n1==NULL)
170     {
171     are1ref=verif_noeud_sur_arete(noeud1->get_coord(),geobloc,&num1,&t1);
172     if (are1ref==NULL)
173     {
174     dim1=2;
175     // face1ref=verif_noeud_sur_face(noeud1->get_coord(),geobloc,&num1);
176     // if (face1ref==NULL) dim1=3; else dim1=2;
177     }
178     else dim1=1;
179     }
180     else
181     {
182     dim1=n1->get_lien_topologie()->get_dimension();
183     if (dim1==1) are1ref=(MG_ARETE*)n1->get_lien_topologie();
184     }
185     if (n2==NULL)
186     {
187     are2ref=verif_noeud_sur_arete(noeud2->get_coord(),geobloc,&num2,&t2);
188     if (are2ref==NULL)
189     {
190     dim2=2;
191     // face2ref=verif_noeud_sur_face(noeud2->get_coord(),geobloc,&num2);
192     // if (face2ref==NULL) dim2=3; else dim2=2;
193     }
194     else dim2=1;
195     }
196     else
197     {
198     dim2=n2->get_lien_topologie()->get_dimension();
199     if (dim2==1) are2ref=(MG_ARETE*)n2->get_lien_topologie();
200     }
201     if ((dim1==3) || (dim2==3)) continue;
202     if ((dim1==2) || (dim2==2))
203     {
204     /*MG_SEGMENT_ADDITIONNEL_FACE segface;
205     segface.segment=seg;
206     segface.noeud1=noeud1;
207     segface.noeud2=noeud2;
208     if (dim1==2) lstsegmentbloc[num1].insert(lstsegmentbloc[num1].end(),segface);
209     else lstsegmentbloc[num2].insert(lstsegmentbloc[num2].end(),segface);*/
210     }
211     else if ((dim1==1) && (dim2==1) && (are1ref!=are2ref))
212     {
213     /*double xyz[3];
214     xyz[0]=0.5*(noeud1->get_coord()[0]+noeud2->get_coord()[0]);
215     xyz[1]=0.5*(noeud1->get_coord()[1]+noeud2->get_coord()[1]);
216     xyz[2]=0.5*(noeud1->get_coord()[2]+noeud2->get_coord()[2]);
217     int num;
218     MG_FACE* face=verif_noeud_sur_face(xyz,geobloc,&num,octreeblocseg);
219     if (face!=NULL)
220     {
221     MG_SEGMENT_ADDITIONNEL_FACE segface;
222     segface.segment=seg;
223     segface.noeud1=noeud1;
224     segface.noeud2=noeud2;
225     lstsegmentbloc[num].insert(lstsegmentbloc[num].end(),segface);
226     } */
227     }
228     else if ((dim1==1) || (dim2==1))
229     {
230     MG_ARETE* ref;
231     if (dim1==1) ref=are1ref; else ref=are2ref;
232     if (n1==NULL)
233     {
234     n1=new MG_NOEUD(ref,noeud1->get_coord()[0],noeud1->get_coord()[1],noeud1->get_coord()[2], IMPOSE);
235     maibloc->ajouter_mg_noeud(n1);
236     octreebloc->inserer(n1);
237     }
238     if (n2==NULL)
239     {
240     n2=new MG_NOEUD(ref,noeud2->get_coord()[0],noeud2->get_coord()[1],noeud2->get_coord()[2], IMPOSE);
241     maibloc->ajouter_mg_noeud(n2);
242     octreebloc->inserer(n2);
243     }
244     MG_SEGMENT *nvsegment=new MG_SEGMENT(ref,n1,n2,IMPOSE);
245     maibloc->ajouter_mg_segment(nvsegment);
246     }
247     else
248     {
249     // cas de 2 noeud sur 2 sommet diff plus tard
250     }
251     }
252     //etape 7 maillage 1D de bloc
253     affiche(" MAILLAGE 1D du modele bloc");
254     int nbaretebloc=geobloc->get_nb_mg_arete();
255     int nbsegmentblocavant=maibloc->get_nb_mg_segment();
256     for (int i=0;i<nbaretebloc;i++)
257     {
258     MAILLEUR1D m1d_bloc(maibloc,geobloc,geobloc->get_mg_arete(i),carte);
259     m1d_bloc.adapte();
260     }
261     // etape 8 transfert des segment vers entier
262     int nbfaceentier=geoentier->get_nb_mg_face();
263     TPL_LISTE_ENTITE<MG_SEGMENT*> *lstsegcontraint=new TPL_LISTE_ENTITE<MG_SEGMENT*>[nbfaceentier];
264     int nbsegmentbloc=maibloc->get_nb_mg_segment();
265     for (int i=nbsegmentblocavant;i<nbsegmentbloc;i++)
266     {
267     MG_SEGMENT* seg=maibloc->get_mg_segment(i);
268     int num;
269     MG_FACE* face=seg_sur_face(seg,geoentier,&num,octreeentierseg);
270     if (face!=NULL)
271     lstsegcontraint[num].ajouter(seg);
272     /*MG_NOEUD* noeud1=seg->get_noeud1();
273     MG_NOEUD* noeud2=seg->get_noeud2();
274     int dim1;
275     int dim2;
276     MG_NOEUD* n1=chercher(noeud1->get_coord(),*octreeentier);
277     MG_NOEUD* n2=chercher(noeud2->get_coord(),*octreeentier);
278     MG_FACE* face1ref;
279     MG_FACE* face2ref;
280     int num1,num2;
281     if (n1==NULL)
282     {
283     face1ref=verif_noeud_sur_face(noeud1->get_coord(),geoentier,&num1,octreeentierseg);
284     if (face1ref==NULL) dim1=3; else dim1=2;
285     }
286     else
287     dim1=n1->get_lien_topologie()->get_dimension();
288     if (n2==NULL)
289     {
290     face2ref=verif_noeud_sur_face(noeud2->get_coord(),geoentier,&num2,octreeentierseg);
291     if (face2ref==NULL) dim2=3; else dim2=2;
292     }
293     else
294     dim2=n2->get_lien_topologie()->get_dimension();
295     if ((dim1==3) || (dim2==3)) continue;
296     if ((dim1==2) || (dim2==2))
297     {
298     int num;
299     if (dim1==2) num=num1;
300     else if (dim2==2) num=num2;
301     lstsegcontraint[num].ajouter(seg);
302     } */
303     }
304     //etape 9 maillage 2D de entier
305     affiche(" MAILLAGE 2D du modele entier");
306     for (int i=0;i<nbfaceentier;i++)
307     {
308     MG_FACE* face=geoentier->get_mg_face(i);
309     MAILLEUR2D m2d(maientier,geoentier,face,carte);
310     m2d.maille(face,&(lstsegcontraint[i]));
311     }
312     //etape 10 transfert des triangle vers bloc
313     int nbfacebloc=geobloc->get_nb_mg_face();
314     TPL_LISTE_ENTITE<MG_TRIANGLE*> *lsttricontraint=new TPL_LISTE_ENTITE<MG_TRIANGLE*>[nbfacebloc];
315     int nbtrientier=maientier->get_nb_mg_triangle();
316     for (int i=0;i<nbtrientier;i++)
317     {
318     MG_TRIANGLE* tri=maientier->get_mg_triangle(i);
319     int num;
320     MG_FACE* face=tri_sur_face(tri,geobloc,&num,octreeblocseg);
321     if (face!=NULL)
322     lsttricontraint[num].ajouter(tri);
323     }
324     //etape 11 maillage 2D de bloc
325     affiche(" MAILLAGE 2D du modele bloc");
326     for (int i=0;i<nbfacebloc;i++)
327     {
328     MG_FACE* face=geobloc->get_mg_face(i);
329     MAILLEUR2D m2d(maibloc,geobloc,face,carte);
330     m2d.maille(face,NULL,&(lsttricontraint[i]));
331     }
332     //etape 12 maillage 3D de bloc
333     affiche(" MAILLAGE 3D du modele bloc");
334     int nbvolumebloc=geobloc->get_nb_mg_volume();
335     for (int i=0;i<nbvolumebloc;i++)
336     {
337     MG_VOLUME *vol=geobloc->get_mg_volume(i);
338     MAILLEUR3D m3d(maibloc,geobloc,vol,carte);
339     m3d.active_affichage(affiche);
340     m3d.maille(vol);
341     }
342     //etape 13 transfert des tetra vers entier
343     TPL_LISTE_ENTITE<MG_TETRA*> lsttetra;
344     LISTE_MG_TETRA::iterator it;
345     for (MG_TETRA* tet=maibloc->get_premier_tetra(it);tet!=NULL;tet=maibloc->get_suivant_tetra(it))
346     lsttetra.ajouter(tet);
347     //etape 14 maillage 3D de entier
348     affiche(" MAILLAGE 3D du modele entier");
349     MG_VOLUME *vol=geoentier->get_mg_volume(0);
350     MAILLEUR3D m3d(maientier,geoentier,vol,carte);
351     m3d.active_affichage(affiche);
352     m3d.maille(vol,NULL,&lsttetra);
353     }
354    
355     MG_NOEUD* MAILLEUR_BLOC::chercher(double *xyz,TPL_OCTREE<MG_NOEUD*,MG_NOEUD*> &octree)
356     {
357     TPL_MAP_ENTITE<MG_NOEUD*> lst;
358     octree.rechercher(xyz[0],xyz[1],xyz[2],1e-5,lst);
359     int nb=lst.get_nb();
360     for (int i=0;i<nb;i++)
361     {
362     OT_VECTEUR_3D vec(xyz,lst.get(i)->get_coord());
363     double metrique[9];
364     carte->evaluer(xyz,metrique);
365     double eps=(1./sqrt(metrique[0])*1e-10);
366     if (vec.get_longueur()<eps) return lst.get(i);
367     }
368     return NULL;
369     }
370    
371    
372     MG_SOMMET* MAILLEUR_BLOC::verif_noeud_sur_sommet(double *xyz, MG_GEOMETRIE *geo)
373     {
374     int nbsom=geo->get_nb_mg_sommet();
375     for (int i=0;i<nbsom;i++)
376     {
377     MG_SOMMET *sommet=geo->get_mg_sommet(i);
378     double xyzsom[3];
379     sommet->get_point()->evaluer(xyzsom);
380     double metrique[9];
381     carte->evaluer(xyz,metrique);
382     double eps=(1./sqrt(metrique[0])*1e-10);
383     OT_VECTEUR_3D vec(xyz,xyzsom);
384     if (vec.get_longueur()<eps) return sommet;
385     }
386     return NULL;
387     }
388    
389    
390     MG_ARETE* MAILLEUR_BLOC::verif_noeud_sur_arete(double *xyz, MG_GEOMETRIE *geo,int *num, double *t)
391     {
392     int nb_are=geo->get_nb_mg_arete();
393     for ( int i=0 ; i < nb_are ; i++)
394     {
395     MG_ARETE *arete=geo->get_mg_arete(i);
396     arete->inverser((*t),xyz);
397     if (arete->get_courbe()->est_periodique())
398     if (*t< arete->get_tmin()) *t=*t+arete->get_courbe()->get_periode();
399     double xyztmp[3];
400     arete->evaluer( *t,xyztmp );
401     double metrique[9];
402     carte->evaluer(xyz,metrique);
403     double eps=(1./sqrt(metrique[0])*1e-10);
404     OT_VECTEUR_3D vec(xyz,xyztmp);
405     if (vec.get_longueur()<eps)
406     if ((*t>arete->get_tmin()) && (*t<arete->get_tmax()))
407     {
408     *num=i;
409     return arete;
410     }
411     }
412     return NULL;
413     }
414    
415     MG_FACE* MAILLEUR_BLOC::seg_sur_face(MG_SEGMENT* seg,MG_GEOMETRIE* geo,int *num,TPL_OCTREE<MG_SEGMENT*,MG_NOEUD*> *octree)
416     {
417     double *xyz1=seg->get_noeud1()->get_coord();
418     double *xyz2=seg->get_noeud2()->get_coord();
419     int nbface=geo->get_nb_mg_face();
420     for (int i=0;i<nbface;i++)
421     {
422     MG_FACE *face=geo->get_mg_face(i);
423     double metrique[9];
424     carte->evaluer(xyz1,metrique);
425     double eps=(1./sqrt(metrique[0])*1e-10);
426     double uv1[2],xyz1bis[3];
427     face->inverser(uv1,xyz1);
428     face->evaluer(uv1,xyz1bis);
429     OT_VECTEUR_3D vec1(xyz1,xyz1bis);
430     if (vec1.get_longueur()<eps)
431     {
432     double uv2[2],xyz2bis[3];
433     face->inverser(uv2,xyz2);
434     face->evaluer(uv2,xyz2bis);
435     OT_VECTEUR_3D vec2(xyz2,xyz2bis);
436     if (vec2.get_longueur()<eps)
437     {
438     double xyz[3];
439     xyz[0]=0.5*(xyz1[0]+xyz2[0]);
440     xyz[1]=0.5*(xyz1[1]+xyz2[1]);
441     xyz[2]=0.5*(xyz1[2]+xyz2[2]);
442     double uv[2];
443     face->inverser(uv,xyz);
444     double dist=calcule_distance_contour_face(xyz,uv,face,octree);
445     if (dist>0.) {*num=i;return face;}
446     }
447     }
448     }
449     return NULL;
450     }
451    
452     MG_FACE* MAILLEUR_BLOC::tri_sur_face(MG_TRIANGLE* tri,MG_GEOMETRIE* geo,int *num,TPL_OCTREE<MG_SEGMENT*,MG_NOEUD*> *octree)
453     {
454     double *xyz1=tri->get_noeud1()->get_coord();
455     double *xyz2=tri->get_noeud2()->get_coord();
456     double *xyz3=tri->get_noeud3()->get_coord();
457     int nbface=geo->get_nb_mg_face();
458     for (int i=0;i<nbface;i++)
459     {
460     MG_FACE *face=geo->get_mg_face(i);
461     double metrique[9];
462     carte->evaluer(xyz1,metrique);
463     double eps=(1./sqrt(metrique[0])*1e-10);
464     double uv1[2],xyz1bis[3];
465     face->inverser(uv1,xyz1);
466     face->evaluer(uv1,xyz1bis);
467     OT_VECTEUR_3D vec1(xyz1,xyz1bis);
468     if (vec1.get_longueur()<eps)
469     {
470     double uv2[2],xyz2bis[3];
471     face->inverser(uv2,xyz2);
472     face->evaluer(uv2,xyz2bis);
473     OT_VECTEUR_3D vec2(xyz2,xyz2bis);
474     if (vec2.get_longueur()<eps)
475     {
476     double uv3[2],xyz3bis[3];
477     face->inverser(uv3,xyz3);
478     face->evaluer(uv3,xyz3bis);
479     OT_VECTEUR_3D vec3(xyz3,xyz3bis);
480     if (vec3.get_longueur()<eps)
481     {
482     double xyz[3];
483     xyz[0]=0.333333333333333333*(xyz1[0]+xyz2[0]+xyz3[0]);
484     xyz[1]=0.333333333333333333*(xyz1[1]+xyz2[1]+xyz3[1]);
485     xyz[2]=0.333333333333333333*(xyz1[2]+xyz2[2]+xyz3[2]);
486     double uv[2];
487     face->inverser(uv,xyz);
488     double dist=calcule_distance_contour_face(xyz,uv,face,octree);
489     if (dist>0.) {*num=i;return face;}
490     }
491     }
492     }
493     }
494     return NULL;
495     }
496    
497     MG_FACE* MAILLEUR_BLOC::verif_noeud_sur_face(double *xyz, MG_GEOMETRIE *geo,int *num,TPL_OCTREE<MG_SEGMENT*,MG_NOEUD*> *octree)
498     {
499     int nbface=geo->get_nb_mg_face();
500     for (int i=0;i<nbface;i++)
501     {
502     MG_FACE *face=geo->get_mg_face(i);
503     double xyz2[3],uv[2];
504     double metrique[9];
505     carte->evaluer(xyz,metrique);
506     double eps=(1./sqrt(metrique[0])*1e-10);
507     face->inverser(uv,xyz);
508     face->evaluer(uv,xyz2);
509     OT_VECTEUR_3D vec(xyz,xyz2);
510     if (vec.get_longueur()<eps)
511     {
512     double dist=calcule_distance_contour_face(xyz,uv,face,octree);
513     if (dist>0.) {*num=i;return face;}
514     }
515     }
516     return NULL;
517     }
518    
519     void MAILLEUR_BLOC::cree_octree(MG_MAILLAGE* mai,TPL_OCTREE<MG_SEGMENT*,MG_NOEUD*> **octreeseg,TPL_OCTREE<MG_NOEUD*,MG_NOEUD*> **octree)
520     {
521     double xmin=1e300,ymin=1e300,zmin=1e300;
522     double xmax=-1e300,ymax=-1e300,zmax=-1e300;
523     TPL_MAP_ENTITE<MG_NOEUD*> lstnoeud;
524     int nb_noeud=mai->get_nb_mg_noeud();
525     for (int i=0;i<nb_noeud;i++)
526     {
527     MG_NOEUD* mgnoeud=mai->get_mg_noeud(i);
528     lstnoeud.ajouter(mgnoeud);
529     xmin=std::min(xmin,mgnoeud->get_x());
530     ymin=std::min(ymin,mgnoeud->get_y());
531     zmin=std::min(zmin,mgnoeud->get_z());
532     xmax=std::max(xmax,mgnoeud->get_x());
533     ymax=std::max(ymax,mgnoeud->get_y());
534     zmax=std::max(zmax,mgnoeud->get_z());
535     }
536     double xdiff=xmax-xmin;
537     double ydiff=ymax-ymin;
538     double zdiff=zmax-zmin;
539     xmin=xmin-0.125*xdiff;
540     xmax=xmax+0.125*xdiff;
541     ymin=ymin-0.125*ydiff;
542     ymax=ymax+0.125*ydiff;
543     zmin=zmin-0.125*zdiff;
544     zmax=zmax+0.125*zdiff;
545     (*octreeseg)=new TPL_OCTREE<MG_SEGMENT*,MG_NOEUD*>;
546     (*octree)=new TPL_OCTREE<MG_NOEUD*,MG_NOEUD*>;
547     (*octreeseg)->initialiser(&lstnoeud,1,xmin,ymin,zmin,xmax,ymax,zmax);
548     (*octree)->initialiser(*octreeseg);
549     int nb_seg=mai->get_nb_mg_segment();
550     for (int i=0;i<nb_seg;i++)
551     (*octreeseg)->inserer(mai->get_mg_segment(i) );
552    
553    
554     }
555    
556    
557    
558    
559     double MAILLEUR_BLOC::calcule_distance_contour_face(double *xyz,double *uv,MG_FACE* face,TPL_OCTREE<MG_SEGMENT*,MG_NOEUD*> *octree)
560     {
561     TPL_MAP_ENTITE<MG_COARETE*> lst;
562     int nbboucle=face->get_nb_mg_boucle();
563     for (int i=0;i<nbboucle;i++)
564     {
565     int nbcoarete=face->get_mg_boucle(i)->get_nb_mg_coarete();
566     for (int j=0;j<nbcoarete;j++)
567     lst.ajouter(face->get_mg_boucle(i)->get_mg_coarete(j));
568     }
569     double metrique[9];
570     carte->evaluer(xyz,metrique);
571     double rayon=1/sqrt(metrique[0]);
572     TPL_MAP_ENTITE<MG_SEGMENT*> lsttrouve;
573     MG_SEGMENT* segretenu=NULL;
574     MG_COARETE* coareteretenu=NULL;
575     double distance=1e300;
576     do
577     {
578     do
579     {
580     octree->rechercher(xyz[0],xyz[1],xyz[2],rayon,lsttrouve);
581     rayon=rayon*1.25;
582     }
583     while (lsttrouve.get_nb()==0);
584     int nb=lsttrouve.get_nb();
585     for (int i=0;i<nb;i++)
586     {
587     MG_SEGMENT* seg=lsttrouve.get(i);
588     MG_ARETE* are=(MG_ARETE*)seg->get_lien_topologie();
589     int nbco=are->get_nb_mg_coarete();
590     MG_COARETE* coarete=NULL;
591     for (int j=0;j<nbco;j++)
592     if (lst.getid(are->get_mg_coarete(j)->get_id())!=NULL) coarete=are->get_mg_coarete(j);
593     if (coarete!=NULL)
594     {
595     double dist=distance_pt_segment(xyz,seg->get_noeud1()->get_coord(),seg->get_noeud2()->get_coord());
596     if (dist<distance)
597     {
598     distance=dist;
599     segretenu=seg;
600     coareteretenu=coarete;
601     }
602     }
603     }
604     }
605     while (segretenu==NULL);
606     OT_VECTEUR_3D base(segretenu->get_noeud1()->get_coord(),segretenu->get_noeud2()->get_coord());
607     base=coareteretenu->get_orientation()*base;
608     double normal[3];
609     face->calcul_normale_unitaire(uv,normal);
610     OT_VECTEUR_3D nor(normal);
611     OT_VECTEUR_3D dir(segretenu->get_noeud1()->get_coord(),xyz);
612     base.norme();
613     dir.norme();
614     double sens=(nor&base)*dir;
615     if (sens<0.) distance=-distance;
616     return distance;
617     }
618    
619    
620     double MAILLEUR_BLOC::distance_pt_segment(double *xyz,double *xyz1,double *xyz2)
621     {
622     double distance;
623     OT_VECTEUR_3D vec12(xyz1,xyz2);
624     OT_VECTEUR_3D vec1p(xyz1,xyz);
625     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());
626     if ((t<0.) || (t>1.))
627     {
628     OT_VECTEUR_3D vec2p(xyz2,xyz);
629     double dist1=vec1p.get_longueur();
630     double dist2=vec2p.get_longueur();
631     distance=0.5*(dist1+dist2);
632     }
633     else
634     {
635     OT_VECTEUR_3D pvec=vec12&vec1p;
636     distance=pvec.get_longueur()/vec12.get_longueur();
637     }
638     return distance;
639     }