ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_bloc.cpp
Revision: 123
Committed: Mon Jul 14 19:13:12 2008 UTC (16 years, 10 months ago) by francois
Original Path: magic/lib/mailleur/mailleur/src/mailleur_bloc.cpp
File size: 26703 byte(s)
Log Message:
calcul des distances dans le parametrique pour le mailleur bloc + ajout d<un vrai cas test

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