ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/optimisation/src/mgopt_posttraitement.cpp
Revision: 224
Committed: Fri Nov 27 15:25:45 2009 UTC (15 years, 5 months ago) by francois
Original Path: magic/lib/optimisation/optimisation/src/mg_lissage.cpp
File size: 16029 byte(s)
Log Message:
Extraction de peau avec detection des cas non manifold pour l'optimisation. Ces cas ne sont pas resolus

File Contents

# User Rev Content
1 francois 222 #include "gestionversion.h"
2     #include "mg_lissage.h"
3     #include "fem_solution.h"
4     #include "fem_maillage.h"
5     #include "mg_maillage.h"
6     #include "mg_gestionnaire.h"
7     #include "mg_triangle_peau.h"
8     #include <stdio.h>
9     #include <string.h>
10 francois 224 #include <math.h>
11 francois 222 //---------------------------------------------------------------------------
12    
13    
14 francois 224
15 francois 222 MG_LISSAGE::MG_LISSAGE()
16     {
17 francois 224 for (int i=0;i<9;i++) etat[i]=0;
18 francois 222 }
19    
20    
21     MG_LISSAGE::~MG_LISSAGE()
22     {
23 francois 224 for (int i=0;i<lst_peau.size();i++)
24     delete lst_peau[i];
25 francois 222 }
26    
27 francois 224 void MG_LISSAGE::conserve(int origine)
28     {
29     etat[(origine-1000)/10]=1;
30     }
31    
32    
33    
34 francois 222 void MG_LISSAGE::extract_skin(FEM_MAILLAGE* mai,MG_GESTIONNAIRE& gest2)
35     {
36    
37     //�tape 0 - On commence par mettre � z�ro tous les nouveau_numero des triangles et des noeuds du maillage
38     MG_MAILLAGE* mg_mai = (MG_MAILLAGE*)mai->get_mg_maillage();
39     LISTE_MG_TRIANGLE::iterator it_tri;
40     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
41     {
42     mgtri->change_nouveau_numero(0);
43 francois 224 mgtri->change_origine(MAILLEUR_AUTO);
44 francois 222 }
45     LISTE_MG_NOEUD::iterator it_noeud;
46     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
47     {
48     mgnoeud->change_nouveau_numero(0);
49     }
50    
51     //�tape 1 - On boucle ensuite tous les t�tra�dres pour ajouter un compteur du nombre de fois qu'un triangle appartient � 1 t�tra
52     LISTE_FEM_TETRA::iterator it_tetra;
53     for (FEM_TETRA* tet=mai->get_premier_tetra(it_tetra);tet!=NULL;tet=mai->get_suivant_tetra(it_tetra))
54     {
55     MG_TETRA* mgtet=(MG_TETRA*)tet->get_mg_element_maillage();
56     int origine = mgtet->get_origine();
57 francois 224 if (origine==IMPOSE)
58 francois 222 {
59 francois 224 mgtet->get_triangle1()->change_origine(IMPOSE);
60     mgtet->get_triangle2()->change_origine(IMPOSE);
61     mgtet->get_triangle3()->change_origine(IMPOSE);
62     mgtet->get_triangle4()->change_origine(IMPOSE);
63     }
64     if (origine==OPTIMISE)
65     {
66     if (mgtet->get_triangle1()->get_origine()!=IMPOSE) mgtet->get_triangle1()->change_origine(OPTIMISE);
67     if (mgtet->get_triangle2()->get_origine()!=IMPOSE) mgtet->get_triangle2()->change_origine(OPTIMISE);
68     if (mgtet->get_triangle3()->get_origine()!=IMPOSE) mgtet->get_triangle3()->change_origine(OPTIMISE);
69     if (mgtet->get_triangle4()->get_origine()!=IMPOSE) mgtet->get_triangle4()->change_origine(OPTIMISE);
70     }
71    
72     if (((origine == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
73     ((origine == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
74     ((origine == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
75     ((origine == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
76     ((origine == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
77     ((origine == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
78    
79     {
80 francois 222 int num1 = mgtet->get_triangle1()->get_nouveau_numero();
81     int num2 = mgtet->get_triangle2()->get_nouveau_numero();
82     int num3 = mgtet->get_triangle3()->get_nouveau_numero();
83     int num4 = mgtet->get_triangle4()->get_nouveau_numero();
84     num1++;
85     num2++;
86     num3++;
87     num4++;
88     mgtet->get_triangle1()->change_nouveau_numero(num1);
89     mgtet->get_triangle2()->change_nouveau_numero(num2);
90     mgtet->get_triangle3()->change_nouveau_numero(num3);
91     mgtet->get_triangle4()->change_nouveau_numero(num4);
92     }
93     }
94    
95     //�tape 2 - On boucle l'ensemble des triangles identifi�s � l'�tape 1 pour identifier les noeuds leur appartenant
96     for (MG_TRIANGLE* mgtri=mg_mai->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mg_mai->get_suivant_triangle(it_tri))
97     {
98     int num = mgtri->get_nouveau_numero();
99     if (num == 1)
100     {
101     MG_NOEUD* noeud1 = mgtri->get_noeud1();
102     MG_NOEUD* noeud2 = mgtri->get_noeud2();
103     MG_NOEUD* noeud3 = mgtri->get_noeud3();
104     noeud1->change_nouveau_numero(1);
105     noeud2->change_nouveau_numero(1);
106     noeud3->change_nouveau_numero(1);
107     }
108     }
109    
110     //�tape 3 - On cr�e un nouveau gestionnaire et un nouveau maillage
111    
112    
113     MG_MAILLAGE* mai2 = new MG_MAILLAGE(NULL);
114     gest2.ajouter_mg_maillage(mai2);
115    
116     // //�tape 4 - On boucle l'ensemble des noeuds identifi�s � l'�tape 2 pour les recr�er dans le second maillage
117     for (MG_NOEUD* mgnoeud=mg_mai->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mg_mai->get_suivant_noeud(it_noeud))
118     {
119     int num = mgnoeud->get_nouveau_numero();
120     if (num == 1)
121     {
122     double x = mgnoeud->get_x();
123     double y = mgnoeud->get_y();
124     double z = mgnoeud->get_z();
125 francois 224 int origine=TRIANGULATION;
126     if (mgnoeud->get_origine()==IMPOSE) origine=IMPOSE;
127     MG_NOEUD* noeud1 = new MG_NOEUD(NULL,x,y,z,origine);
128 francois 222 mai2->ajouter_mg_noeud(noeud1);
129     mgnoeud->change_nouveau_numero(noeud1->get_id());
130     }
131     }
132    
133     // //�tape 5 - On boucle l'ensemble des triangles identifi�s � l'�tape 1 pour les recr�er dans le maillage 2
134 francois 224 for (FEM_TETRA* tet=mai->get_premier_tetra(it_tetra);tet!=NULL;tet=mai->get_suivant_tetra(it_tetra))
135 francois 222 {
136 francois 224 MG_TETRA* mgtet=(MG_TETRA*)tet->get_mg_element_maillage();
137     int originetet=mgtet->get_origine();
138     if (((originetet == OPTIMISE) && (etat[(OPTIMISE-1000)/10]==1) ) ||
139     ((originetet == IMPOSE) && (etat[(IMPOSE-1000)/10]==1) ) ||
140     ((originetet == MAILLEUR_AUTO) && (etat[(MAILLEUR_AUTO-1000)/10]==1) ) ||
141     ((originetet == TRIANGULATION) && (etat[(TRIANGULATION-1000)/10]==1) ) ||
142     ((originetet == MODIFICATION) && (etat[(MODIFICATION-1000)/10]==1) ) ||
143     ((originetet == DUPLIQUER) && (etat[(DUPLIQUER-1000)/10]==1) ) )
144 francois 222 {
145 francois 224 MG_NOEUD* noeud1 = mgtet->get_noeud1();
146     MG_NOEUD* noeud2 = mgtet->get_noeud2();
147     MG_NOEUD* noeud3 = mgtet->get_noeud3();
148     MG_NOEUD* noeud4 = mgtet->get_noeud4();
149 francois 222 MG_NOEUD* node1 = mai2->get_mg_noeudid(noeud1->get_nouveau_numero());
150     MG_NOEUD* node2 = mai2->get_mg_noeudid(noeud2->get_nouveau_numero());
151     MG_NOEUD* node3 = mai2->get_mg_noeudid(noeud3->get_nouveau_numero());
152 francois 224 MG_NOEUD* node4 = mai2->get_mg_noeudid(noeud4->get_nouveau_numero());
153     MG_TRIANGLE* tri1=mgtet->get_triangle1();
154     MG_TRIANGLE* tri2=mgtet->get_triangle2();
155     MG_TRIANGLE* tri3=mgtet->get_triangle3();
156     MG_TRIANGLE* tri4=mgtet->get_triangle4();
157     if (tri1->get_nouveau_numero()==1)
158     {
159     int origine=TRIANGULATION;
160     if (tri1->get_origine()==IMPOSE) origine=IMPOSE;
161     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node2,node3,mai2,origine);
162     tripeau->change_nouveau_numero(0);
163     tri1->change_nouveau_numero(0);
164     }
165     if (tri2->get_nouveau_numero()==1)
166     {
167     int origine=TRIANGULATION;
168     if (tri2->get_origine()==IMPOSE) origine=IMPOSE;
169     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node4,node2,mai2,origine);
170     tripeau->change_nouveau_numero(0);
171     tri2->change_nouveau_numero(0);
172     }
173     if (tri3->get_nouveau_numero()==1)
174     {
175     int origine=TRIANGULATION;
176     if (tri3->get_origine()==IMPOSE) origine=IMPOSE;
177     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node2,node4,node3,mai2,origine);
178     tripeau->change_nouveau_numero(0);
179     tri3->change_nouveau_numero(0);
180     }
181     if (tri4->get_nouveau_numero()==1)
182     {
183     int origine=TRIANGULATION;
184     if (tri4->get_origine()==IMPOSE) origine=IMPOSE;
185     MG_TRIANGLE_PEAU* tripeau = insere_triangle(NULL,node1,node3,node4,mai2,origine);
186     tripeau->change_nouveau_numero(0);
187     tri4->change_nouveau_numero(0);
188     }
189 francois 222 }
190 francois 224 }
191     cout << mai2->get_nb_mg_triangle()<< " triangles" <<endl;
192     // etape 6 recherche des vosins
193     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
194     {
195     MG_TRIANGLE_PEAU* tripeau=(MG_TRIANGLE_PEAU*)mgtri;
196     MG_TRIANGLE_PEAU* voisin1=recherche_voisin(tripeau->get_noeud1(),tripeau->get_noeud2(),tripeau);
197     MG_TRIANGLE_PEAU* voisin2=recherche_voisin(tripeau->get_noeud2(),tripeau->get_noeud3(),tripeau);
198     MG_TRIANGLE_PEAU* voisin3=recherche_voisin(tripeau->get_noeud3(),tripeau->get_noeud1(),tripeau);
199     tripeau->change_voisin1(voisin1);
200     tripeau->change_voisin2(voisin2);
201     tripeau->change_voisin3(voisin3);
202     }
203     // //�tape 7 - On recherche les peaux
204     int fin;
205     do
206     {
207     fin=1;
208     for (MG_TRIANGLE* mgtri=mai2->get_premier_triangle(it_tri);mgtri!=NULL;mgtri=mai2->get_suivant_triangle(it_tri))
209     {
210     MG_TRIANGLE_PEAU *tripeau=(MG_TRIANGLE_PEAU*)mgtri;
211     if (tripeau->get_nouveau_numero()==0)
212     {
213     fin=0;
214     std::vector<MG_TRIANGLE_PEAU*> *peau=new std::vector<MG_TRIANGLE_PEAU*>;
215     lst_peau.push_back(peau);
216     tripeau->change_nouveau_numero(1);
217     peau->push_back(tripeau);
218     determine_peau(peau);
219     }
220     }
221     }
222     while (fin==0);
223     std::cout << lst_peau.size() << " peau" << endl;
224     for (int i=0;i<lst_peau.size();i++)
225     std::cout << " peau " << i << " " << lst_peau[i]->size() << " triangles" << endl;
226     //test de manifold
227     for (MG_NOEUD* mgnoeud=mai2->get_premier_noeud(it_noeud);mgnoeud!=NULL;mgnoeud=mai2->get_suivant_noeud(it_noeud))
228     mgnoeud->change_nouveau_numero(0);
229     LISTE_MG_SEGMENT::iterator itseg;
230     for (MG_SEGMENT* seg=mai2->get_premier_segment(itseg);seg!=NULL;seg=mai2->get_suivant_segment(itseg))
231     seg->change_nouveau_numero(0);
232     TPL_MAP_ENTITE<MG_SEGMENT*> nonmanifoldarete;
233     TPL_MAP_ENTITE<MG_NOEUD*> nonmanifoldnoeud;
234     int nbpeau=lst_peau.size();
235     for (int i=0;i<nbpeau;i++)
236     {
237     int nbtri=lst_peau[i]->size();
238     for (int j=0;j<nbtri;j++)
239     {
240     MG_TRIANGLE_PEAU* tripeau=(*lst_peau[i])[j];
241     tripeau->get_segment1()->change_nouveau_numero(tripeau->get_segment1()->get_nouveau_numero()+1);
242     tripeau->get_segment2()->change_nouveau_numero(tripeau->get_segment2()->get_nouveau_numero()+1);
243     tripeau->get_segment3()->change_nouveau_numero(tripeau->get_segment3()->get_nouveau_numero()+1);
244     if (tripeau->get_segment1()->get_nouveau_numero()>2)
245     nonmanifoldarete.ajouter(tripeau->get_segment1());
246     if (tripeau->get_segment2()->get_nouveau_numero()>2)
247     nonmanifoldarete.ajouter(tripeau->get_segment2());
248     if (tripeau->get_segment3()->get_nouveau_numero()>2)
249     nonmanifoldarete.ajouter(tripeau->get_segment3());
250     if ((tripeau->get_noeud1()->get_nouveau_numero()!=0) && (tripeau->get_noeud1()->get_nouveau_numero()!=i+1))
251     nonmanifoldnoeud.ajouter(tripeau->get_noeud1());
252     if ((tripeau->get_noeud2()->get_nouveau_numero()!=0) && (tripeau->get_noeud2()->get_nouveau_numero()!=i+1))
253     nonmanifoldnoeud.ajouter(tripeau->get_noeud2());
254     if ((tripeau->get_noeud3()->get_nouveau_numero()!=0) && (tripeau->get_noeud3()->get_nouveau_numero()!=i+1))
255     nonmanifoldnoeud.ajouter(tripeau->get_noeud3());
256     tripeau->get_noeud1()->change_nouveau_numero(i+1);
257     tripeau->get_noeud2()->change_nouveau_numero(i+1);
258     tripeau->get_noeud3()->change_nouveau_numero(i+1);
259     }
260     }
261     int nbnonmanifoldarete=nonmanifoldarete.get_nb();
262     for (int i=0;i<nbnonmanifoldarete;i++)
263     {
264     MG_SEGMENT* seg=nonmanifoldarete.get(i);
265     MG_NOEUD* n1=seg->get_noeud1();
266     MG_NOEUD* n2=seg->get_noeud2();
267     nonmanifoldnoeud.supprimer(n1);
268     nonmanifoldnoeud.supprimer(n2);
269     }
270 francois 222
271 francois 224
272     cout << "Non manifold par arete " << nonmanifoldarete.get_nb() << endl;
273     cout << "Non manifold par noeud " << nonmanifoldnoeud.get_nb() << endl;;
274 francois 222 }
275    
276 francois 224
277     void MG_LISSAGE::determine_peau(std::vector<MG_TRIANGLE_PEAU*> * peau)
278 francois 222 {
279 francois 224 int num=0;
280     while (num!=peau->size())
281     {
282     MG_TRIANGLE_PEAU* p=(*peau)[num];
283     if (p->get_voisin1()->get_nouveau_numero()==0)
284     {
285     peau->push_back(p->get_voisin1());
286     p->get_voisin1()->change_nouveau_numero(1);
287     }
288     if (p->get_voisin2()->get_nouveau_numero()==0)
289     {
290     peau->push_back(p->get_voisin2());
291     p->get_voisin2()->change_nouveau_numero(1);
292     }
293     if (p->get_voisin3()->get_nouveau_numero()==0)
294     {
295     peau->push_back(p->get_voisin3());
296     p->get_voisin3()->change_nouveau_numero(1);
297     }
298     num++;
299     }
300     }
301    
302    
303     MG_TRIANGLE_PEAU* MG_LISSAGE::insere_triangle(MG_ELEMENT_TOPOLOGIQUE* topo,class MG_NOEUD *mgnoeud1,class MG_NOEUD *mgnoeud2,class MG_NOEUD *mgnoeud3,MG_MAILLAGE* mg_maillage,int origine)
304     {
305 francois 222 MG_TRIANGLE_PEAU* voisin1=NULL,*voisin2=NULL,*voisin3=NULL;
306     MG_SEGMENT* mgsegment1=mg_maillage->get_mg_segment(mgnoeud1->get_id(),mgnoeud2->get_id());
307     MG_SEGMENT* mgsegment2=mg_maillage->get_mg_segment(mgnoeud2->get_id(),mgnoeud3->get_id());
308     MG_SEGMENT* mgsegment3=mg_maillage->get_mg_segment(mgnoeud3->get_id(),mgnoeud1->get_id());
309     if (mgsegment1==NULL)
310 francois 224 mgsegment1=mg_maillage->ajouter_mg_segment(topo,mgnoeud1,mgnoeud2,origine);
311 francois 222 if (mgsegment2==NULL)
312 francois 224 mgsegment2=mg_maillage->ajouter_mg_segment(topo,mgnoeud2,mgnoeud3,origine);
313 francois 222 if (mgsegment3==NULL)
314 francois 224 mgsegment3=mg_maillage->ajouter_mg_segment(topo,mgnoeud3,mgnoeud1,origine);
315     MG_TRIANGLE_PEAU* mtriangle=new MG_TRIANGLE_PEAU(topo,mgnoeud1,mgnoeud2,mgnoeud3,mgsegment1,mgsegment2,mgsegment3,origine);
316 francois 222 mg_maillage->ajouter_mg_triangle(mtriangle);
317     return mtriangle;
318     }
319    
320 francois 224 MG_TRIANGLE_PEAU* MG_LISSAGE::recherche_voisin(MG_NOEUD* mg_noeud1,MG_NOEUD* mg_noeud2,MG_TRIANGLE_PEAU* triref)
321 francois 222 {
322 francois 224 MG_TRIANGLE_PEAU* trisol=NULL;
323     double angleref=4.*M_PI;
324 francois 222 int nb1=mg_noeud1->get_lien_triangle()->get_nb();
325     int nb2=mg_noeud2->get_lien_triangle()->get_nb();
326     for (int i=0;i<nb1;i++)
327     for (int j=0;j<nb2;j++)
328 francois 224 if (mg_noeud1->get_lien_triangle()->get(i)==mg_noeud2->get_lien_triangle()->get(j))
329     if ((MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i)!=triref)
330     {
331     double angle=calcul_angle(triref,(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i));
332     if (angle<angleref)
333     {
334     angleref=angle;
335     trisol=(MG_TRIANGLE_PEAU*)mg_noeud1->get_lien_triangle()->get(i);
336     }
337     }
338     return trisol;
339 francois 222 }
340 francois 223
341 francois 224
342    
343     double MG_LISSAGE::calcul_angle(MG_TRIANGLE_PEAU* ft1,MG_TRIANGLE_PEAU* ft2)
344 francois 223 {
345 francois 224 MG_NOEUD* noeud1=ft1->get_noeud1();
346     MG_NOEUD* noeud2=ft1->get_noeud2();
347     MG_NOEUD* noeud3=ft1->get_noeud3();
348     MG_NOEUD* noeud4=ft2->get_noeud1();
349     MG_NOEUD* noeud5=ft2->get_noeud2();
350     MG_NOEUD* noeud6=ft2->get_noeud3();
351    
352     MG_NOEUD* noeuda=NULL;
353     MG_NOEUD* noeudb=NULL;
354     MG_NOEUD* noeudc=NULL;
355     MG_NOEUD* noeudd=NULL;
356    
357    
358     if (!(noeud1==noeud4))
359     if (!(noeud1==noeud5))
360     if (!(noeud1==noeud6))
361     {
362     noeuda=noeud1;
363     noeudb=noeud2;
364     noeudc=noeud3;
365     }
366     if (!(noeud2==noeud4))
367     if (!(noeud2==noeud5))
368     if (!(noeud2==noeud6))
369     {
370     noeuda=noeud2;
371     noeudb=noeud1;
372     noeudc=noeud3;
373     }
374     if (!(noeud3==noeud4))
375     if (!(noeud3==noeud5))
376     if (!(noeud3==noeud6))
377     {
378     noeuda=noeud3;
379     noeudb=noeud1;
380     noeudc=noeud2;
381     }
382     if (noeuda==NULL) return 2*M_PI;
383     if (!(noeud4==noeudb))
384     if (!(noeud4==noeudc)) noeudd=noeud4;
385     if (!(noeud5==noeudb))
386     if (!(noeud5==noeudc)) noeudd=noeud5;
387     if (!(noeud6==noeudb))
388     if (!(noeud6==noeudc)) noeudd=noeud6;
389    
390     OT_VECTEUR_3D milieu(0.5*(noeudb->get_x()+noeudc->get_x()),0.5*(noeudb->get_y()+noeudc->get_y()),0.5*(noeudb->get_z()+noeudc->get_z()));
391     OT_VECTEUR_3D bc(noeudc->get_x()-noeudb->get_x(),noeudc->get_y()-noeudb->get_y(),noeudc->get_z()-noeudb->get_z());
392     bc.norme();
393     double xyz[3];
394     xyz[0]=noeuda->get_x();
395     xyz[1]=noeuda->get_y();
396     xyz[2]=noeuda->get_z();
397     OT_VECTEUR_3D n1(xyz,milieu.get_xyz());
398     n1.norme();
399     OT_VECTEUR_3D n=bc&n1;
400     n1=n&bc;
401     n1.norme();
402     xyz[0]=noeudd->get_x();
403     xyz[1]=noeudd->get_y();
404     xyz[2]=noeudd->get_z();
405     OT_VECTEUR_3D n2(xyz,milieu.get_xyz());
406     n2.norme();
407     n=bc&n2;
408     n2=n&bc;
409     n2.norme();
410    
411     double ps=n1*n2;
412     if (ps>1.) ps=1.;
413     if (ps<-1.) ps=-1.;
414     double angle=acos(ps);
415     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
416     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
417     OT_VECTEUR_3D nbnd(noeudd->get_x()-noeudb->get_x(),noeudd->get_y()-noeudb->get_y(),noeudd->get_z()-noeudb->get_z());
418     n=n1n2&n1n3;
419     n.norme();
420     nbnd.norme();
421    
422     ps=n*nbnd;
423     if (ps<-0.0001) angle=2*M_PI-angle;
424     return angle;
425    
426    
427 francois 223 }