ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/remailleur.cpp
Revision: 165
Committed: Wed Oct 8 21:45:56 2008 UTC (16 years, 7 months ago) by francois
Original Path: magic/lib/mailleur/mailleur/src/remailleur.cpp
File size: 50899 byte(s)
Log Message:
Resolution D'un bug sur le pb d'un point appartient a une face. Je suis pas sur que ca marche tout le temps. Mais jusqu au prochain bug c est bon

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