ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/mg_arete.cpp
Revision: 576
Committed: Wed Oct 22 18:13:01 2014 UTC (10 years, 6 months ago) by francois
File size: 14494 byte(s)
Log Message: 
ajout d'une geometrie virtuelle dans le maillage structure pour pouvoir saisir des conditions aux limites. ajout de conditions aux limites pour le thermique

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 francois 576 bool MG_ARETE::est_une_arete_element(void)
110     {
111     return false;
112     }
113 francois 283
114     void MG_ARETE::evaluer(double t,double *xyz)
115     {
116     if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
117     courbe->evaluer(t,xyz);
118     }
119    
120     void MG_ARETE::deriver(double t,double *xyz)
121     {
122     if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
123     courbe->deriver(t,xyz);
124     if (orientation!=MEME_SENS)
125     {
126     xyz[0]=-xyz[0];
127     xyz[1]=-xyz[1];
128     xyz[2]=-xyz[2];
129     }
130     }
131    
132     void MG_ARETE::deriver_seconde(double t,double *ddxyz,double* dxyz,double* xyz)
133     {
134     if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
135     courbe->deriver_seconde(t,ddxyz,dxyz,xyz);
136     if (orientation!=MEME_SENS)
137     {
138     dxyz[0]=-dxyz[0];
139     dxyz[1]=-dxyz[1];
140     dxyz[2]=-dxyz[2];
141     }
142     }
143    
144     void MG_ARETE::inverser(double& t,double *xyz,double precision)
145     {
146     courbe->inverser(t,xyz,precision);
147     if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
148     }
149    
150     double MG_ARETE::get_tmin(void)
151     {
152     return cosommet1->get_t();
153     }
154    
155     double MG_ARETE::get_tmax(void)
156     {
157     if (!courbe->est_periodique()) return cosommet2->get_t();
158     double tmin=cosommet1->get_t();
159     double tmax=cosommet2->get_t();
160     if (tmax<tmin+1e-6) tmax=tmax+courbe->get_periode();
161     return tmax;
162     }
163    
164     double MG_ARETE::get_longueur(double t1,double t2,double precis)
165     {
166     if (orientation!=MEME_SENS)
167     {
168     t1=courbe->get_tmin()+courbe->get_tmax()-t1;
169     t2=courbe->get_tmin()+courbe->get_tmax()-t2;
170     return (-courbe->get_longueur(t1,t2,precis));
171     }
172     return courbe->get_longueur(t1,t2,precis);
173    
174     }
175    
176     double MG_ARETE::get_M(double t)
177     {
178     if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
179     return courbe->get_M(t);
180     }
181    
182    
183    
184     void MG_ARETE::supprimer_mg_coarete(class MG_COARETE* coarete)
185     {
186     std::vector<MG_COARETE*>::iterator i;
187     for (i=lst_coarete.begin();i!=lst_coarete.end();i++)
188     {
189     if ((*i)==coarete)
190     {
191     lst_coarete.erase(i);
192     return;
193     }
194     }
195     }
196    
197    
198    
199     MG_COARETE* MG_ARETE::get_mg_coarete(int num)
200     {
201     return lst_coarete[num];
202     }
203    
204    
205 gervaislavoie 332 TPL_LISTE_ENTITE<class MG_SEGMENT*>* MG_ARETE::get_lien_segment(void)
206 francois 283 {
207     return &segment;
208     }
209    
210     int MG_ARETE::get_dimension(void)
211     {
212     return 1;
213     }
214    
215    
216     VCT& MG_ARETE::get_vectorisation(void)
217     {
218     if (vect==NULL) vect=new VCT_ARETE(this);
219     return *vect;
220     }
221    
222     void MG_ARETE::enregistrer(std::ostream& o)
223     {
224     int nb=get_nb_ccf();
225     o << "%" << get_id() << "=ARETE("<< get_idoriginal() << ",$" << courbe->get_id() << ",$"<<cosommet1->get_id() << ",$" <<cosommet2->get_id() << "," << orientation << "," << nb;
226     if (nb!=0)
227     {
228     o << ",(";
229     for (int i=0;i<nb;i++)
230     {
231     char nom[3];
232     get_type_ccf(i,nom);
233     o << "(" << nom << "," << get_valeur_ccf(i) << ")";
234     if (i!=nb-1) o << "," ;
235     }
236     o << ")";
237     }
238     o << ");" << std::endl;
239     }
240    
241    
242     void MG_ARETE::get_intersection(double t1, double t2,double* point_iners)
243     {
244     OT_VECTEUR_3D P1_INTERSECTION;
245    
246     double U,W;
247    
248     double xyz2[3];
249     double xyz1[3];
250     double drv1[3];
251     double drv2[3];
252    
253     evaluer(t1,xyz1);
254     deriver(t1,drv1);
255    
256     evaluer(t2,xyz2);
257     deriver(t2,drv2);
258    
259     OT_VECTEUR_3D a(xyz1) ;
260     OT_VECTEUR_3D b(drv1) ;
261     OT_VECTEUR_3D c(xyz2) ;
262     OT_VECTEUR_3D d(drv2) ;
263     OT_VECTEUR_3D f=c-a;
264     d=-d;
265     double det1=b[0]*d[1]-b[1]*d[0];
266     double det2=b[0]*d[2]-b[2]*d[0];
267     double det3=b[1]*d[2]-b[2]*d[1];
268    
269     bool c_fait=false;
270    
271     if (fabs(det1)>=1e-10)
272     {
273     double u,w;
274     u=(d[1]*f[0]-d[0]*f[1])/det1;
275     w=(b[0]*f[1]-b[1]*f[0])/det1;
276     double verif_x=a[2] +u*b[2];
277     double verif_y=c[2]-w*d[2];
278     if (fabs(verif_x-verif_y)<0.000001)
279     {
280     U=u ;
281     W=w ;
282     c_fait=true;
283     }
284    
285     }
286    
287     if (!c_fait&&fabs(det2)>=1e-10)
288     {
289     double u,w;
290     u=(d[2]*f[0]-d[0]*f[2])/det2;
291     w=(b[0]*f[2]-b[2]*f[0])/det2;
292     double verif_x=a[1]+u*b[1];
293     double verif_y=c[1]-w*d[1];
294     if (fabs(verif_x-verif_y)<0.000001)
295     {
296     U=u ;
297     W=w ;
298     c_fait=true;
299     }
300     }
301    
302     if (!c_fait&&fabs(det3)>=1e-10)
303     {
304     double u,w;
305     u=(d[2]*f[1]-d[1]*f[2])/det3;
306     w=(b[1]*f[2]-b[2]*f[1])/det3;
307     double verif_x=a[0] +u*b[0];
308     double verif_y=c[0]-w*d[0];
309     if (fabs(verif_x-verif_y)<0.000001)
310     {
311     U=u ;
312     W=w ;
313     c_fait=true;
314     }
315     }
316    
317     P1_INTERSECTION=a+U*b;
318    
319     point_iners[0]=P1_INTERSECTION[0];
320     point_iners[1]=P1_INTERSECTION[1];
321     point_iners[2]=P1_INTERSECTION[2];
322    
323     }
324    
325    
326    
327     void MG_ARETE:: get_param_NURBS(int& indx_premier_ptctr,TPL_LISTE_ENTITE<double> &param)
328     {
329    
330     TPL_LISTE_ENTITE<double> param1;
331     int type=courbe->get_type_geometrique(param1);
332    
333     double tcs1=get_tmin();
334     double tcs2=get_tmax();
335     double angle=tcs1-tcs2;
336    
337     if ((type==MGCo_CIRCLE||type==MGCo_ELLIPSE)&&(tcs1!=tcs2)&&fabs(angle)!=2*PI)
338     {
339    
340     double xyz[3];
341     param.ajouter(1);
342     param.ajouter(4);
343     param.ajouter(0);
344     param.ajouter(7);
345     param.ajouter(0);
346    
347     double arc_lengh= get_longueur(tcs1,tcs2);
348     double rayon=arc_lengh/fabs(angle);
349    
350     double dif=fabs(fabs(angle)-PI);
351     if (fabs(fabs(angle)-PI)<=1e-6||fabs(angle)-PI<0.) //ouverture inf ou egale a PI
352     {
353    
354    
355     double request_lengh=arc_lengh/2;
356    
357     param.ajouter(0);
358     param.ajouter(0);
359     param.ajouter(0);
360     param.ajouter(0.5);
361     param.ajouter(0.5);
362     param.ajouter(1);
363     param.ajouter(1);
364     param.ajouter(1);
365    
366    
367     double xyz1[3];
368     double xyz1_inters[3];
369     double xyz2_inters[3];
370     double xyz2[3];
371     double xyz3[3];
372     double t_mid=get_t(tcs1,tcs2,request_lengh);
373    
374     evaluer(tcs1,xyz1);
375    
376     param.ajouter(xyz1[0]);
377     param.ajouter(xyz1[1]);
378     param.ajouter(xyz1[2]);
379     param.ajouter(1);
380    
381     get_intersection(tcs1, t_mid,xyz1_inters);
382    
383     param.ajouter(xyz1_inters[0]);
384     param.ajouter(xyz1_inters[1]);
385     param.ajouter(xyz1_inters[2]);
386     param.ajouter(0.5);
387    
388     get_intersection(t_mid, tcs2,xyz2_inters);
389    
390     param.ajouter(xyz2_inters[0]);
391     param.ajouter(xyz2_inters[1]);
392     param.ajouter(xyz2_inters[2]);
393     param.ajouter(0.5);
394    
395     evaluer(tcs2,xyz3);
396    
397     param.ajouter(xyz3[0]);
398     param.ajouter(xyz3[1]);
399     param.ajouter(xyz3[2]);
400     param.ajouter(1);
401    
402     indx_premier_ptctr=13;
403    
404     }
405    
406     if (fabs(fabs(angle)-PI)>1e-6&&fabs(angle)-PI>0.)
407     {
408    
409     double xyz[3];
410    
411 francois 290 double lengh_at_mid=arc_lengh;
412     double request_lengh=lengh_at_mid/2.;
413     double tmid=get_t(tcs1,tcs2,request_lengh);
414     double tmid1=get_t(tcs1,tmid,request_lengh/2.);
415     double tmid2=get_t(tmid,tcs2,request_lengh/2.);
416 francois 283
417     param.ajouter(0);
418     param.ajouter(0);
419     param.ajouter(0);
420     param.ajouter(0.25);
421     param.ajouter(0.25);
422     param.ajouter(0.5);
423     param.ajouter(0.5);
424     param.ajouter(0.75);
425     param.ajouter(0.75);
426     param.ajouter(1);
427     param.ajouter(1);
428     param.ajouter(1);
429    
430    
431    
432     evaluer(tcs1,xyz);
433    
434     param.ajouter(xyz[0]);
435     param.ajouter(xyz[1]);
436     param.ajouter(xyz[2]);
437     param.ajouter(1);
438    
439 francois 290
440     get_intersection(tcs1,tmid1,xyz);
441 francois 283
442     param.ajouter(xyz[0]);
443     param.ajouter(xyz[1]);
444     param.ajouter(xyz[2]);
445     param.ajouter(0.5);
446    
447 francois 290 get_intersection(tmid1, tmid,xyz);
448 francois 283
449     param.ajouter(xyz[0]);
450     param.ajouter(xyz[1]);
451     param.ajouter(xyz[2]);
452     param.ajouter(0.5);
453    
454 francois 290 evaluer(tmid,xyz);
455 francois 283
456     param.ajouter(xyz[0]);
457     param.ajouter(xyz[1]);
458     param.ajouter(xyz[2]);
459     param.ajouter(1);
460    
461 francois 290
462     get_intersection(tmid,tmid2,xyz);
463 francois 283
464     param.ajouter(xyz[0]);
465     param.ajouter(xyz[1]);
466     param.ajouter(xyz[2]);
467     param.ajouter(0.5);
468    
469    
470 francois 290 get_intersection(tmid2,tcs2,xyz);
471 francois 283
472     param.ajouter(xyz[0]);
473     param.ajouter(xyz[1]);
474     param.ajouter(xyz[2]);
475     param.ajouter(0.5);
476    
477     evaluer(tcs2,xyz);
478    
479     param.ajouter(xyz[0]);
480     param.ajouter(xyz[1]);
481     param.ajouter(xyz[2]);
482     param.ajouter(1);
483    
484     indx_premier_ptctr=17;
485     }
486    
487     }
488    
489    
490     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    
493     if (type==MGCo_BSPLINE)
494     courbe->get_param_NURBS(indx_premier_ptctr,param);
495    
496     if (type==MGCo_LINE)
497     {
498    
499     param.ajouter(1);
500     param.ajouter(2);
501     param.ajouter(0);
502     param.ajouter(2);
503     param.ajouter(0);
504    
505    
506     param.ajouter(0);
507     param.ajouter(0);
508     param.ajouter(1);
509     param.ajouter(1);
510    
511     double xyz1[3];
512     double xyz2[3];
513    
514     double tcs1=get_tmin();
515     double tcs2=get_tmax();
516    
517     evaluer(tcs1,xyz1);
518    
519     param.ajouter(xyz1[0]);
520     param.ajouter(xyz1[1]);
521     param.ajouter(xyz1[2]);
522     param.ajouter(1);
523    
524     evaluer(tcs2,xyz2);
525    
526     param.ajouter(xyz2[0]);
527     param.ajouter(xyz2[1]);
528     param.ajouter(xyz2[2]);
529     param.ajouter(1);
530    
531     indx_premier_ptctr=9;
532     }
533    
534     }
535    
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     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     }
584