ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/mg_arete.cpp
Revision: 77
Committed: Wed Apr 2 23:55:19 2008 UTC (17 years, 1 month ago) by francois
Original Path: magic/lib/geometrie/geometrie/src/mg_arete.cpp
File size: 12682 byte(s)
Log Message: 
interpolation lineaire dans le get milieu + precision ajuster dans double2

File Contents

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