ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/mg_arete.cpp
Revision: 375
Committed: Fri Nov 9 15:27:14 2012 UTC (12 years, 6 months ago) by francois
File size: 14433 byte(s)
Log Message: 
Solution FEM sont vectorielles et tensorielles + affichage des options de compilations dans les exe

File Contents

# User Rev Content
1 francois 283 //------------------------------------------------------------
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     #include "vct_arete.h"
28     //#include "message.h"
29     //#include "affiche.h"
30 francois 375 #include "mg_definition.h"
31 francois 283 #include <math.h>
32     #include "ot_mathematique.h"
33     #include "constantegeo.h"
34     #define PI 3.1415926535897932384626433832795
35    
36     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     {
38     }
39    
40     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     {
42     }
43    
44     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     {
46     }
47    
48     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     {
50     }
51     MG_ARETE::MG_ARETE(std::string idori,class MG_COURBE* crb,int sens):MG_ELEMENT_TOPOLOGIQUE(idori),courbe(crb),orientation(sens),vect(NULL)
52     {
53     }
54     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    
64     MG_ARETE::~MG_ARETE()
65     {
66     // if (lst_coarete.size()!=0) afficheur << WARCOARETE << this->get_id()<< enderr;
67     segment.vide();
68     if (vect!=NULL) delete vect;
69     }
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 gervaislavoie 332 TPL_LISTE_ENTITE<class MG_SEGMENT*>* MG_ARETE::get_lien_segment(void)
202 francois 283 {
203     return &segment;
204     }
205    
206     int MG_ARETE::get_dimension(void)
207     {
208     return 1;
209     }
210    
211    
212     VCT& MG_ARETE::get_vectorisation(void)
213     {
214     if (vect==NULL) vect=new VCT_ARETE(this);
215     return *vect;
216     }
217    
218     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     void MG_ARETE::get_intersection(double t1, double t2,double* point_iners)
239     {
240     OT_VECTEUR_3D P1_INTERSECTION;
241    
242     double U,W;
243    
244     double xyz2[3];
245     double xyz1[3];
246     double drv1[3];
247     double drv2[3];
248    
249     evaluer(t1,xyz1);
250     deriver(t1,drv1);
251    
252     evaluer(t2,xyz2);
253     deriver(t2,drv2);
254    
255     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    
265     bool c_fait=false;
266    
267     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    
281     }
282    
283     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    
298     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    
313     P1_INTERSECTION=a+U*b;
314    
315     point_iners[0]=P1_INTERSECTION[0];
316     point_iners[1]=P1_INTERSECTION[1];
317     point_iners[2]=P1_INTERSECTION[2];
318    
319     }
320    
321    
322    
323     void MG_ARETE:: get_param_NURBS(int& indx_premier_ptctr,TPL_LISTE_ENTITE<double> &param)
324     {
325    
326     TPL_LISTE_ENTITE<double> param1;
327     int type=courbe->get_type_geometrique(param1);
328    
329     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     {
335    
336     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     double arc_lengh= get_longueur(tcs1,tcs2);
344     double rayon=arc_lengh/fabs(angle);
345    
346     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    
350    
351     double request_lengh=arc_lengh/2;
352    
353     param.ajouter(0);
354     param.ajouter(0);
355     param.ajouter(0);
356     param.ajouter(0.5);
357     param.ajouter(0.5);
358     param.ajouter(1);
359     param.ajouter(1);
360     param.ajouter(1);
361    
362    
363     double xyz1[3];
364     double xyz1_inters[3];
365     double xyz2_inters[3];
366     double xyz2[3];
367     double xyz3[3];
368     double t_mid=get_t(tcs1,tcs2,request_lengh);
369    
370     evaluer(tcs1,xyz1);
371    
372     param.ajouter(xyz1[0]);
373     param.ajouter(xyz1[1]);
374     param.ajouter(xyz1[2]);
375     param.ajouter(1);
376    
377     get_intersection(tcs1, t_mid,xyz1_inters);
378    
379     param.ajouter(xyz1_inters[0]);
380     param.ajouter(xyz1_inters[1]);
381     param.ajouter(xyz1_inters[2]);
382     param.ajouter(0.5);
383    
384     get_intersection(t_mid, tcs2,xyz2_inters);
385    
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     evaluer(tcs2,xyz3);
392    
393     param.ajouter(xyz3[0]);
394     param.ajouter(xyz3[1]);
395     param.ajouter(xyz3[2]);
396     param.ajouter(1);
397    
398     indx_premier_ptctr=13;
399    
400     }
401    
402     if (fabs(fabs(angle)-PI)>1e-6&&fabs(angle)-PI>0.)
403     {
404    
405     double xyz[3];
406    
407 francois 290 double lengh_at_mid=arc_lengh;
408     double request_lengh=lengh_at_mid/2.;
409     double tmid=get_t(tcs1,tcs2,request_lengh);
410     double tmid1=get_t(tcs1,tmid,request_lengh/2.);
411     double tmid2=get_t(tmid,tcs2,request_lengh/2.);
412 francois 283
413     param.ajouter(0);
414     param.ajouter(0);
415     param.ajouter(0);
416     param.ajouter(0.25);
417     param.ajouter(0.25);
418     param.ajouter(0.5);
419     param.ajouter(0.5);
420     param.ajouter(0.75);
421     param.ajouter(0.75);
422     param.ajouter(1);
423     param.ajouter(1);
424     param.ajouter(1);
425    
426    
427    
428     evaluer(tcs1,xyz);
429    
430     param.ajouter(xyz[0]);
431     param.ajouter(xyz[1]);
432     param.ajouter(xyz[2]);
433     param.ajouter(1);
434    
435 francois 290
436     get_intersection(tcs1,tmid1,xyz);
437 francois 283
438     param.ajouter(xyz[0]);
439     param.ajouter(xyz[1]);
440     param.ajouter(xyz[2]);
441     param.ajouter(0.5);
442    
443 francois 290 get_intersection(tmid1, tmid,xyz);
444 francois 283
445     param.ajouter(xyz[0]);
446     param.ajouter(xyz[1]);
447     param.ajouter(xyz[2]);
448     param.ajouter(0.5);
449    
450 francois 290 evaluer(tmid,xyz);
451 francois 283
452     param.ajouter(xyz[0]);
453     param.ajouter(xyz[1]);
454     param.ajouter(xyz[2]);
455     param.ajouter(1);
456    
457 francois 290
458     get_intersection(tmid,tmid2,xyz);
459 francois 283
460     param.ajouter(xyz[0]);
461     param.ajouter(xyz[1]);
462     param.ajouter(xyz[2]);
463     param.ajouter(0.5);
464    
465    
466 francois 290 get_intersection(tmid2,tcs2,xyz);
467 francois 283
468     param.ajouter(xyz[0]);
469     param.ajouter(xyz[1]);
470     param.ajouter(xyz[2]);
471     param.ajouter(0.5);
472    
473     evaluer(tcs2,xyz);
474    
475     param.ajouter(xyz[0]);
476     param.ajouter(xyz[1]);
477     param.ajouter(xyz[2]);
478     param.ajouter(1);
479    
480     indx_premier_ptctr=17;
481     }
482    
483     }
484    
485    
486     if ((type==MGCo_ELLIPSE||type==MGCo_CIRCLE)&&OPERATEUR::egal(fabs(angle),2*PI,0.000001))
487     courbe->get_param_NURBS(indx_premier_ptctr,param);
488    
489     if (type==MGCo_BSPLINE)
490     courbe->get_param_NURBS(indx_premier_ptctr,param);
491    
492     if (type==MGCo_LINE)
493     {
494    
495     param.ajouter(1);
496     param.ajouter(2);
497     param.ajouter(0);
498     param.ajouter(2);
499     param.ajouter(0);
500    
501    
502     param.ajouter(0);
503     param.ajouter(0);
504     param.ajouter(1);
505     param.ajouter(1);
506    
507     double xyz1[3];
508     double xyz2[3];
509    
510     double tcs1=get_tmin();
511     double tcs2=get_tmax();
512    
513     evaluer(tcs1,xyz1);
514    
515     param.ajouter(xyz1[0]);
516     param.ajouter(xyz1[1]);
517     param.ajouter(xyz1[2]);
518     param.ajouter(1);
519    
520     evaluer(tcs2,xyz2);
521    
522     param.ajouter(xyz2[0]);
523     param.ajouter(xyz2[1]);
524     param.ajouter(xyz2[2]);
525     param.ajouter(1);
526    
527     indx_premier_ptctr=9;
528     }
529    
530     }
531    
532     struct integrale {
533     double ti;
534     double li;
535     };
536    
537     double MG_ARETE:: get_t(double t1,double t2,double longueur_voulue,double dt)
538     {
539     double pt[3];
540     double t,ti,tii;
541    
542     //double longueur_voulue,longueur=0;
543     double longueur=0;
544     std::vector<integrale> tab_integrale;
545     //longueur_voulue= alpha *get_longueur(t1,t2);
546    
547     integrale pas;
548     ti=t1;
549     do
550     {
551     tii=ti+(t2-t1)*dt;
552    
553     t=0.7886751345*ti+0.2113248654*tii;
554     deriver(t,pt);
555    
556     double facteur=pow(pt[0],2)+pow(pt[1],2)+pow(pt[2],2);
557     longueur=longueur+0.5*(tii-ti)*sqrt(facteur);
558     t=0.7886751345*tii+0.2113248654*ti;
559     deriver(t,pt);
560     facteur=pow(pt[0],2)+pow(pt[1],2)+pow(pt[2],2);
561     longueur=longueur+0.5*(tii-ti)*sqrt(facteur);
562     integrale pas;
563     pas.ti=ti;
564     pas.li=longueur;
565     tab_integrale.insert(tab_integrale.end(),pas);
566     ti=tii;
567     }
568     while (fabs(tab_integrale.back().li)<fabs(longueur_voulue));
569    
570     //double t_trouve=0.7886751345*tab_integrale.back().ti+0.2113248654*tab_integrale.back().ti+(t2-t1)*dt;
571     double tn_1=tab_integrale.back().ti;
572     double ln_1=tab_integrale[tab_integrale.size()-2].li;
573     double tn=ti;
574     double ln=tab_integrale.back().li;
575    
576    
577     double t_trouve=tn_1+(tn-tn_1)/(ln-ln_1)*(longueur_voulue-ln_1);
578     return t_trouve ;
579     }
580