ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/remailleur.cpp
Revision: 168
Committed: Wed Dec 10 19:53:59 2008 UTC (16 years, 5 months ago) by francois
Original Path: magic/lib/mailleur/mailleur/src/remailleur.cpp
File size: 51652 byte(s)
Log Message:
Nouvelle zone de destruction pour le remailleur. Optimisation meilleure.

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