ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur_bloc.cpp
Revision: 129
Committed: Thu Jul 24 20:44:37 2008 UTC (16 years, 9 months ago) by francois
Original Path: magic/lib/mailleur/mailleur/src/mailleur_bloc.cpp
File size: 32698 byte(s)
Log Message:
octree mieux raffinee dans le mailleur 3D + accent ote dans le texte pour meilleur compatibilité avec linux 

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