ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/mg_arete.cpp
Revision: 253
Committed: Tue Jul 13 19:40:46 2010 UTC (14 years, 10 months ago) by francois
File size: 13008 byte(s)
Log Message: 
changement de hiearchie et utilisation de ccmake + mise a jour

File Contents

# User Rev Content
1 5 //------------------------------------------------------------
2     //------------------------------------------------------------
3     // MAGiC
4 francois 195 // Jean Christophe Cuilli�re et Vincent FRANCOIS
5     // D�partement de G�nie M�canique - UQTR
6 5 //------------------------------------------------------------
7 francois 195 // Le projet MAGIC est un projet de recherche du d�partement
8     // de g�nie m�canique de l'Universit� du Qu�bec �
9     // Trois Rivi�res
10     // Les librairies ne peuvent �tre utilis�es sans l'accord
11 5 // des auteurs (contact : francois@uqtr.ca)
12     //------------------------------------------------------------
13     //------------------------------------------------------------
14     //
15     // mg_arete.cpp
16     //
17     //------------------------------------------------------------
18     //------------------------------------------------------------
19     // COPYRIGHT 2000
20 francois 195 // Version du 02/03/2006 � 11H22
21 5 //------------------------------------------------------------
22     //------------------------------------------------------------
23    
24    
25     #include "gestionversion.h"
26     #include "mg_arete.h"
27 souaissa 71 #include "vct_arete.h"
28 5 //#include "message.h"
29     //#include "affiche.h"
30     #include "geom.h"
31     #include <math.h>
32     #include "ot_mathematique.h"
33     #include "constantegeo.h"
34 souaissa 73 #define PI 3.1415926535897932384626433832795
35 5
36 souaissa 71 MG_ARETE::MG_ARETE(std::string idori,unsigned long num,MG_COSOMMET* mgcosom1,MG_COSOMMET* mgcosom2,MG_COURBE* crb,int sens):MG_ELEMENT_TOPOLOGIQUE(num,idori),cosommet1(mgcosom1),cosommet2(mgcosom2),courbe(crb),orientation(sens),vect(NULL)
37 5 {
38     }
39    
40 souaissa 71 MG_ARETE::MG_ARETE(std::string idori,MG_COSOMMET* mgcosom1,MG_COSOMMET* mgcosom2,MG_COURBE* crb,int sens):MG_ELEMENT_TOPOLOGIQUE(idori),cosommet1(mgcosom1),cosommet2(mgcosom2),courbe(crb),orientation(sens),vect(NULL)
41 5 {
42     }
43    
44 souaissa 71 MG_ARETE::MG_ARETE(MG_ARETE& mdd):MG_ELEMENT_TOPOLOGIQUE(mdd),cosommet1(mdd.cosommet1),cosommet2(mdd.cosommet2),courbe(mdd.courbe),orientation(mdd.orientation),vect(NULL)
45 5 {
46     }
47    
48 souaissa 71 MG_ARETE::MG_ARETE(std::string idori,unsigned long num,class MG_COURBE* crb,int sens):MG_ELEMENT_TOPOLOGIQUE(num,idori),courbe(crb),orientation(sens),vect(NULL)
49 5 {
50     }
51 souaissa 71 MG_ARETE::MG_ARETE(std::string idori,class MG_COURBE* crb,int sens):MG_ELEMENT_TOPOLOGIQUE(idori),courbe(crb),orientation(sens),vect(NULL)
52 5 {
53     }
54 francois 195 void MG_ARETE::get_topologie_sousjacente(TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> *lst)
55     {
56     MG_SOMMET* som1=cosommet1->get_sommet();
57     MG_SOMMET* som2=cosommet2->get_sommet();
58     lst->ajouter(som1);
59     lst->ajouter(som2);
60     som1->get_topologie_sousjacente(lst);
61     som2->get_topologie_sousjacente(lst);
62     }
63 5
64     MG_ARETE::~MG_ARETE()
65     {
66     // if (lst_coarete.size()!=0) afficheur << WARCOARETE << this->get_id()<< enderr;
67     segment.vide();
68 souaissa 71 if (vect!=NULL) delete vect;
69 5 }
70    
71     void MG_ARETE::changer_cosommet1(class MG_COSOMMET* cosom)
72     {
73     cosommet1=cosom;
74     }
75     void MG_ARETE::changer_cosommet2(class MG_COSOMMET* cosom)
76     {
77     cosommet2=cosom;
78     }
79     MG_COSOMMET* MG_ARETE::get_cosommet1(void)
80     {
81     return cosommet1;
82     }
83     MG_COSOMMET* MG_ARETE::get_cosommet2(void)
84     {
85     return cosommet2;
86     }
87     MG_COURBE* MG_ARETE::get_courbe(void)
88     {
89     return courbe;
90     }
91     int MG_ARETE::get_orientation(void)
92     {
93     return orientation;
94     }
95    
96    
97    
98    
99     void MG_ARETE::ajouter_mg_coarete(class MG_COARETE* coarete)
100     {
101     lst_coarete.insert(lst_coarete.end(),coarete);
102     }
103    
104     int MG_ARETE::get_nb_mg_coarete(void)
105     {
106     return lst_coarete.size();
107     }
108    
109    
110     void MG_ARETE::evaluer(double t,double *xyz)
111     {
112     if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
113     courbe->evaluer(t,xyz);
114     }
115    
116     void MG_ARETE::deriver(double t,double *xyz)
117     {
118     if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
119     courbe->deriver(t,xyz);
120     if (orientation!=MEME_SENS)
121     {
122     xyz[0]=-xyz[0];
123     xyz[1]=-xyz[1];
124     xyz[2]=-xyz[2];
125     }
126     }
127    
128     void MG_ARETE::deriver_seconde(double t,double *ddxyz,double* dxyz,double* xyz)
129     {
130     if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
131     courbe->deriver_seconde(t,ddxyz,dxyz,xyz);
132     if (orientation!=MEME_SENS)
133     {
134     dxyz[0]=-dxyz[0];
135     dxyz[1]=-dxyz[1];
136     dxyz[2]=-dxyz[2];
137     }
138     }
139    
140     void MG_ARETE::inverser(double& t,double *xyz,double precision)
141     {
142     courbe->inverser(t,xyz,precision);
143     if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
144     }
145    
146     double MG_ARETE::get_tmin(void)
147     {
148     return cosommet1->get_t();
149     }
150    
151     double MG_ARETE::get_tmax(void)
152     {
153     if (!courbe->est_periodique()) return cosommet2->get_t();
154     double tmin=cosommet1->get_t();
155     double tmax=cosommet2->get_t();
156     if (tmax<tmin+1e-6) tmax=tmax+courbe->get_periode();
157     return tmax;
158     }
159    
160     double MG_ARETE::get_longueur(double t1,double t2,double precis)
161     {
162     if (orientation!=MEME_SENS)
163     {
164     t1=courbe->get_tmin()+courbe->get_tmax()-t1;
165     t2=courbe->get_tmin()+courbe->get_tmax()-t2;
166     return (-courbe->get_longueur(t1,t2,precis));
167     }
168     return courbe->get_longueur(t1,t2,precis);
169    
170     }
171    
172     double MG_ARETE::get_M(double t)
173     {
174     if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
175     return courbe->get_M(t);
176     }
177    
178    
179    
180     void MG_ARETE::supprimer_mg_coarete(class MG_COARETE* coarete)
181     {
182     std::vector<MG_COARETE*>::iterator i;
183     for (i=lst_coarete.begin();i!=lst_coarete.end();i++)
184     {
185     if ((*i)==coarete)
186     {
187     lst_coarete.erase(i);
188     return;
189     }
190     }
191     }
192    
193    
194    
195     MG_COARETE* MG_ARETE::get_mg_coarete(int num)
196     {
197     return lst_coarete[num];
198     }
199    
200    
201     TPL_LISTE_ENTITE<class MG_SEGMENT*>* MG_ARETE::get_lien_segement(void)
202     {
203     return &segment;
204     }
205    
206     int MG_ARETE::get_dimension(void)
207     {
208     return 1;
209     }
210    
211 souaissa 71
212     VCT& MG_ARETE::get_vectorisation(void)
213     {
214     if (vect==NULL) vect=new VCT_ARETE(this);
215     return *vect;
216     }
217    
218 5 void MG_ARETE::enregistrer(std::ostream& o)
219     {
220     int nb=get_nb_ccf();
221     o << "%" << get_id() << "=ARETE("<< get_idoriginal() << ",$" << courbe->get_id() << ",$"<<cosommet1->get_id() << ",$" <<cosommet2->get_id() << "," << orientation << "," << nb;
222     if (nb!=0)
223     {
224     o << ",(";
225     for (int i=0;i<nb;i++)
226     {
227     char nom[3];
228     get_type_ccf(i,nom);
229     o << "(" << nom << "," << get_valeur_ccf(i) << ")";
230     if (i!=nb-1) o << "," ;
231     }
232     o << ")";
233     }
234     o << ");" << std::endl;
235     }
236    
237    
238 souaissa 73 void MG_ARETE::get_intersection(double t1, double t2,double* point_iners)
239 5 {
240 souaissa 73 OT_VECTEUR_3D P1_INTERSECTION;
241    
242     double U,W;
243    
244     double xyz2[3];
245 5 double xyz1[3];
246 souaissa 73 double drv1[3];
247     double drv2[3];
248 5
249 francois 19 evaluer(t1,xyz1);
250 souaissa 73 deriver(t1,drv1);
251 5
252 francois 19 evaluer(t2,xyz2);
253 souaissa 73 deriver(t2,drv2);
254 5
255 souaissa 73 OT_VECTEUR_3D a(xyz1) ;
256     OT_VECTEUR_3D b(drv1) ;
257     OT_VECTEUR_3D c(xyz2) ;
258     OT_VECTEUR_3D d(drv2) ;
259     OT_VECTEUR_3D f=c-a;
260     d=-d;
261     double det1=b[0]*d[1]-b[1]*d[0];
262     double det2=b[0]*d[2]-b[2]*d[0];
263     double det3=b[1]*d[2]-b[2]*d[1];
264 5
265 souaissa 73 bool c_fait=false;
266 5
267 souaissa 73 if(fabs(det1)>=1e-10)
268     {
269     double u,w;
270     u=(d[1]*f[0]-d[0]*f[1])/det1;
271     w=(b[0]*f[1]-b[1]*f[0])/det1;
272     double verif_x=a[2] +u*b[2];
273     double verif_y=c[2]-w*d[2];
274     if(fabs(verif_x-verif_y)<0.000001)
275     {
276     U=u ;
277     W=w ;
278     c_fait=true;
279     }
280 5
281 souaissa 73 }
282 5
283 souaissa 73 if (!c_fait&&fabs(det2)>=1e-10)
284     {
285     double u,w;
286     u=(d[2]*f[0]-d[0]*f[2])/det2;
287     w=(b[0]*f[2]-b[2]*f[0])/det2;
288     double verif_x=a[1]+u*b[1];
289     double verif_y=c[1]-w*d[1];
290     if(fabs(verif_x-verif_y)<0.000001)
291     {
292     U=u ;
293     W=w ;
294     c_fait=true;
295     }
296     }
297 5
298 souaissa 73 if (!c_fait&&fabs(det3)>=1e-10)
299     {
300     double u,w;
301     u=(d[2]*f[1]-d[1]*f[2])/det3;
302     w=(b[1]*f[2]-b[2]*f[1])/det3;
303     double verif_x=a[0] +u*b[0];
304     double verif_y=c[0]-w*d[0];
305     if(fabs(verif_x-verif_y)<0.000001)
306     {
307     U=u ;
308     W=w ;
309     c_fait=true;
310     }
311     }
312 5
313 souaissa 73 P1_INTERSECTION=a+U*b;
314 5
315 francois 19 point_iners[0]=P1_INTERSECTION[0];
316     point_iners[1]=P1_INTERSECTION[1];
317     point_iners[2]=P1_INTERSECTION[2];
318 5
319 francois 19 }
320 5
321    
322    
323 francois 19 void MG_ARETE:: get_param_NURBS(int& indx_premier_ptctr,TPL_LISTE_ENTITE<double> &param)
324     {
325 5
326 francois 19 TPL_LISTE_ENTITE<double> param1;
327     int type=courbe->get_type_geometrique(param1);
328    
329 souaissa 73 double tcs1=get_tmin();
330     double tcs2=get_tmax();
331     double angle=tcs1-tcs2;
332    
333     if((type==MGCo_CIRCLE||type==MGCo_ELLIPSE)&&(tcs1!=tcs2)&&fabs(angle)!=2*PI)
334 francois 19 {
335    
336 5 double xyz[3];
337     param.ajouter(1);
338     param.ajouter(4);
339     param.ajouter(0);
340     param.ajouter(7);
341     param.ajouter(0);
342    
343 souaissa 73 double arc_lengh= get_longueur(tcs1,tcs2);
344     double rayon=arc_lengh/fabs(angle);
345 5
346 souaissa 73 double dif=fabs(fabs(angle)-PI);
347     if (fabs(fabs(angle)-PI)<=1e-6||fabs(angle)-PI<0.) //ouverture inf ou egale a PI
348     {
349 francois 19
350    
351 souaissa 73 double request_lengh=arc_lengh/2;
352 francois 19
353 5 param.ajouter(0);
354     param.ajouter(0);
355     param.ajouter(0);
356 francois 19 param.ajouter(0.5);
357     param.ajouter(0.5);
358 5 param.ajouter(1);
359     param.ajouter(1);
360     param.ajouter(1);
361 francois 19
362    
363     double xyz1[3];
364     double xyz1_inters[3];
365 souaissa 73 double xyz2_inters[3];
366     double xyz2[3];
367 francois 19 double xyz3[3];
368 souaissa 73 double t_mid=get_t(tcs1,tcs2,request_lengh);
369 francois 19
370     evaluer(tcs1,xyz1);
371    
372     param.ajouter(xyz1[0]);
373     param.ajouter(xyz1[1]);
374     param.ajouter(xyz1[2]);
375 5 param.ajouter(1);
376    
377 souaissa 73 get_intersection(tcs1, t_mid,xyz1_inters);
378 5
379 francois 19 param.ajouter(xyz1_inters[0]);
380     param.ajouter(xyz1_inters[1]);
381     param.ajouter(xyz1_inters[2]);
382     param.ajouter(0.5);
383    
384 souaissa 73 get_intersection(t_mid, tcs2,xyz2_inters);
385 francois 19
386     param.ajouter(xyz2_inters[0]);
387     param.ajouter(xyz2_inters[1]);
388     param.ajouter(xyz2_inters[2]);
389     param.ajouter(0.5);
390    
391 souaissa 73 evaluer(tcs2,xyz3);
392 francois 19
393     param.ajouter(xyz3[0]);
394     param.ajouter(xyz3[1]);
395 5 param.ajouter(xyz3[2]);
396     param.ajouter(1);
397    
398 francois 19 indx_premier_ptctr=13;
399    
400     }
401 souaissa 73
402     if (fabs(fabs(angle)-PI)>1e-6&&fabs(angle)-PI>0.)
403 francois 19 {
404    
405 souaissa 73 double xyz[3];
406 francois 19
407 souaissa 73 double lengh_at_pi=rayon*PI;
408     double request_lengh=0.5*rayon*PI;
409     double t1_milieu=get_t(tcs1,tcs2,request_lengh);
410    
411 francois 19 param.ajouter(0);
412     param.ajouter(0);
413     param.ajouter(0);
414     param.ajouter(0.25);
415     param.ajouter(0.25);
416 5 param.ajouter(0.5);
417 francois 19 param.ajouter(0.5);
418     param.ajouter(0.75);
419     param.ajouter(0.75);
420     param.ajouter(1);
421     param.ajouter(1);
422     param.ajouter(1);
423 5
424    
425 francois 19
426 souaissa 73 evaluer(tcs1,xyz);
427    
428     param.ajouter(xyz[0]);
429     param.ajouter(xyz[1]);
430     param.ajouter(xyz[2]);
431 francois 19 param.ajouter(1);
432    
433 souaissa 73 double t2_milieu=get_t(t1_milieu,tcs2,request_lengh);
434 francois 19
435 souaissa 73 get_intersection(tcs1,t1_milieu,xyz);
436 francois 19
437 souaissa 73 param.ajouter(xyz[0]);
438     param.ajouter(xyz[1]);
439     param.ajouter(xyz[2]);
440 5 param.ajouter(0.5);
441    
442 souaissa 73 get_intersection(t1_milieu, t2_milieu,xyz);
443    
444     param.ajouter(xyz[0]);
445     param.ajouter(xyz[1]);
446     param.ajouter(xyz[2]);
447 francois 19 param.ajouter(0.5);
448    
449 souaissa 73 evaluer(t2_milieu,xyz);
450 francois 19
451 souaissa 73 param.ajouter(xyz[0]);
452     param.ajouter(xyz[1]);
453     param.ajouter(xyz[2]);
454 5 param.ajouter(1);
455 francois 19
456 souaissa 73 request_lengh=arc_lengh-lengh_at_pi;
457     request_lengh=0.5*request_lengh;
458 francois 19
459 souaissa 73 t1_milieu=get_t(t2_milieu,tcs2,request_lengh);
460 francois 19
461    
462 souaissa 73 get_intersection(t2_milieu,t1_milieu,xyz);
463 francois 19
464 souaissa 73 param.ajouter(xyz[0]);
465     param.ajouter(xyz[1]);
466     param.ajouter(xyz[2]);
467     param.ajouter(0.5);
468 francois 19
469    
470 souaissa 73 get_intersection(t1_milieu,tcs2,xyz);
471 francois 19
472 souaissa 73 param.ajouter(xyz[0]);
473     param.ajouter(xyz[1]);
474     param.ajouter(xyz[2]);
475 francois 19 param.ajouter(0.5);
476    
477 souaissa 73 evaluer(tcs2,xyz);
478 francois 19
479 souaissa 73 param.ajouter(xyz[0]);
480     param.ajouter(xyz[1]);
481     param.ajouter(xyz[2]);
482 francois 19 param.ajouter(1);
483    
484     indx_premier_ptctr=17;
485     }
486 souaissa 73
487 5 }
488    
489    
490 souaissa 73 if ((type==MGCo_ELLIPSE||type==MGCo_CIRCLE)&&OPERATEUR::egal(fabs(angle),2*PI,0.000001))
491     courbe->get_param_NURBS(indx_premier_ptctr,param);
492 5
493 souaissa 73 if (type==MGCo_BSPLINE)
494 francois 19 courbe->get_param_NURBS(indx_premier_ptctr,param);
495 5
496 souaissa 73 if(type==MGCo_LINE)
497 5 {
498    
499     param.ajouter(1);
500     param.ajouter(2);
501 francois 19 param.ajouter(0);
502 5 param.ajouter(2);
503 francois 19 param.ajouter(0);
504 5
505    
506     param.ajouter(0);
507     param.ajouter(0);
508     param.ajouter(1);
509     param.ajouter(1);
510    
511 souaissa 73 double xyz1[3];
512     double xyz2[3];
513 5
514 souaissa 73 double tcs1=get_tmin();
515     double tcs2=get_tmax();
516 5
517 souaissa 73 evaluer(tcs1,xyz1);
518 5
519 souaissa 73 param.ajouter(xyz1[0]);
520     param.ajouter(xyz1[1]);
521     param.ajouter(xyz1[2]);
522 5 param.ajouter(1);
523    
524 souaissa 73 evaluer(tcs2,xyz2);
525 5
526 souaissa 73 param.ajouter(xyz2[0]);
527     param.ajouter(xyz2[1]);
528     param.ajouter(xyz2[2]);
529 5 param.ajouter(1);
530 francois 19
531     indx_premier_ptctr=9;
532 5 }
533    
534     }
535 souaissa 73
536     struct integrale{
537     double ti;
538     double li;
539     };
540    
541     double MG_ARETE:: get_t(double t1,double t2,double longueur_voulue,double dt)
542     {
543     double pt[3];
544     double t,ti,tii;
545    
546     //double longueur_voulue,longueur=0;
547     double longueur=0;
548     std::vector<integrale> tab_integrale;
549     //longueur_voulue= alpha *get_longueur(t1,t2);
550    
551     integrale pas;
552     ti=t1;
553     do
554     {
555     tii=ti+(t2-t1)*dt;
556    
557     t=0.7886751345*ti+0.2113248654*tii;
558     deriver(t,pt);
559    
560     double facteur=pow(pt[0],2)+pow(pt[1],2)+pow(pt[2],2);
561     longueur=longueur+0.5*(tii-ti)*sqrt(facteur);
562     t=0.7886751345*tii+0.2113248654*ti;
563     deriver(t,pt);
564     facteur=pow(pt[0],2)+pow(pt[1],2)+pow(pt[2],2);
565     longueur=longueur+0.5*(tii-ti)*sqrt(facteur);
566     integrale pas;
567     pas.ti=ti;
568     pas.li=longueur;
569     tab_integrale.insert(tab_integrale.end(),pas);
570     ti=tii;
571     }
572     while(fabs(tab_integrale.back().li)<fabs(longueur_voulue));
573    
574     //double t_trouve=0.7886751345*tab_integrale.back().ti+0.2113248654*tab_integrale.back().ti+(t2-t1)*dt;
575 francois 77 double tn_1=tab_integrale.back().ti;
576     double ln_1=tab_integrale[tab_integrale.size()-2].li;
577     double tn=ti;
578     double ln=tab_integrale.back().li;
579    
580    
581     double t_trouve=tn_1+(tn-tn_1)/(ln-ln_1)*(longueur_voulue-ln_1);
582     return t_trouve ;
583 souaissa 73 }
584