ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/mg_arete.cpp
Revision: 906
Committed: Mon Nov 13 22:30:18 2017 UTC (7 years, 9 months ago) by couturad
File size: 17595 byte(s)
Log Message:
Nouveau opencascade commit 1

File Contents

# Content
1 //------------------------------------------------------------
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 #include "mg_definition.h"
31 #include <math.h>
32 #include "ot_mathematique.h"
33 #include "constantegeo.h"
34 #include "fem_solution.h"
35 #include "fem_maillage.h"
36 #define PI 3.1415926535897932384626433832795
37
38 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)
39 {
40 }
41
42 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)
43 {
44 }
45
46 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)
47 {
48 }
49
50 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)
51 {
52 cosommet1=NULL;
53 cosommet2=NULL;
54 }
55 MG_ARETE::MG_ARETE(std::string idori,class MG_COURBE* crb,int sens):MG_ELEMENT_TOPOLOGIQUE(idori),courbe(crb),orientation(sens),vect(NULL)
56 {
57 cosommet1=NULL;
58 cosommet2=NULL;
59 }
60 void MG_ARETE::get_topologie_sousjacente(TPL_MAP_ENTITE<MG_ELEMENT_TOPOLOGIQUE*> *lst)
61 {
62 MG_SOMMET* som1=cosommet1->get_sommet();
63 MG_SOMMET* som2=cosommet2->get_sommet();
64 lst->ajouter(som1);
65 lst->ajouter(som2);
66 som1->get_topologie_sousjacente(lst);
67 som2->get_topologie_sousjacente(lst);
68 }
69
70 MG_ARETE::~MG_ARETE()
71 {
72 // if (lst_coarete.size()!=0) afficheur << WARCOARETE << this->get_id()<< enderr;
73 if (vect!=NULL) delete vect;
74 }
75
76 void MG_ARETE::changer_cosommet1(class MG_COSOMMET* cosom)
77 {
78 cosommet1=cosom;
79 }
80 void MG_ARETE::changer_cosommet2(class MG_COSOMMET* cosom)
81 {
82 cosommet2=cosom;
83 }
84 MG_COSOMMET* MG_ARETE::get_cosommet1(void)
85 {
86 return cosommet1;
87 }
88 MG_COSOMMET* MG_ARETE::get_cosommet2(void)
89 {
90 return cosommet2;
91 }
92 MG_COURBE* MG_ARETE::get_courbe(void)
93 {
94 return courbe;
95 }
96 int MG_ARETE::get_orientation(void)
97 {
98 return orientation;
99 }
100
101
102
103
104 void MG_ARETE::ajouter_mg_coarete(class MG_COARETE* coarete)
105 {
106 lst_coarete.insert(lst_coarete.end(),coarete);
107 }
108
109 int MG_ARETE::get_nb_mg_coarete(void)
110 {
111 return lst_coarete.size();
112 }
113
114 bool MG_ARETE::est_une_arete_element(void)
115 {
116 return false;
117 }
118
119 void MG_ARETE::evaluer(double t,double *xyz)
120 {
121 if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
122 courbe->evaluer(t,xyz);
123 }
124
125 void MG_ARETE::deriver(double t,double *xyz)
126 {
127 if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
128 courbe->deriver(t,xyz);
129 if (orientation!=MEME_SENS)
130 {
131 xyz[0]=-xyz[0];
132 xyz[1]=-xyz[1];
133 xyz[2]=-xyz[2];
134 }
135 }
136
137 void MG_ARETE::deriver_seconde(double t,double *ddxyz,double* dxyz,double* xyz)
138 {
139 if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
140 courbe->deriver_seconde(t,ddxyz,dxyz,xyz);
141 if (orientation!=MEME_SENS)
142 {
143 dxyz[0]=-dxyz[0];
144 dxyz[1]=-dxyz[1];
145 dxyz[2]=-dxyz[2];
146 }
147 }
148
149 void MG_ARETE::inverser(double& t,double *xyz,double precision)
150 {
151 courbe->inverser(t,xyz,precision);
152 if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
153 }
154
155 double MG_ARETE::get_tmin(void)
156 {
157 return cosommet1->get_t();
158 }
159
160 double MG_ARETE::get_tmax(void)
161 {
162 if (!courbe->est_periodique()) return cosommet2->get_t();
163 double tmin=cosommet1->get_t();
164 double tmax=cosommet2->get_t();
165 if (tmax<tmin+1e-6) tmax=tmax+courbe->get_periode();
166 return tmax;
167 }
168
169 double MG_ARETE::get_longueur(double t1,double t2,double precis)
170 {
171 if (orientation!=MEME_SENS)
172 {
173 t1=courbe->get_tmin()+courbe->get_tmax()-t1;
174 t2=courbe->get_tmin()+courbe->get_tmax()-t2;
175 return (-courbe->get_longueur(t1,t2,precis));
176 }
177 return courbe->get_longueur(t1,t2,precis);
178
179 }
180
181 double MG_ARETE::get_M(double t)
182 {
183 if (orientation!=MEME_SENS) t=courbe->get_tmin()+courbe->get_tmax()-t;
184 return courbe->get_M(t);
185 }
186
187
188
189 void MG_ARETE::supprimer_mg_coarete(class MG_COARETE* coarete)
190 {
191 std::vector<MG_COARETE*>::iterator i;
192 for (i=lst_coarete.begin();i!=lst_coarete.end();i++)
193 {
194 if ((*i)==coarete)
195 {
196 lst_coarete.erase(i);
197 return;
198 }
199 }
200 }
201
202
203
204 MG_COARETE* MG_ARETE::get_mg_coarete(int num)
205 {
206 return lst_coarete[num];
207 }
208
209 int MG_ARETE::get_type(void)
210 {
211 return TYPE_ELEMENT_TOPOLOGIQUE::ARETE;
212 }
213
214
215 int MG_ARETE::get_dimension(void)
216 {
217 return 1;
218 }
219
220
221 VCT& MG_ARETE::get_vectorisation(void)
222 {
223 if (vect==NULL) vect=new VCT_ARETE(this);
224 return *vect;
225 }
226
227 void MG_ARETE::enregistrer(std::ostream& o,double version)
228 {
229 int nb=get_nb_ccf();
230 o << "%" << get_id() << "=ARETE("<< get_idoriginal() << ",$" << courbe->get_id() << ",$"<<cosommet1->get_id() << ",$" <<cosommet2->get_id() << "," << orientation << "," ;
231 if (version<2)
232 {
233 o<< nb;
234 if (nb!=0)
235 {
236 o << ",(";
237 for (int i=0;i<nb;i++)
238 {
239 char nom[3];
240 get_type_ccf(i,nom);
241 o << "(" << nom << "," << get_valeur_ccf(i) << ")";
242 if (i!=nb-1) o << "," ;
243 }
244 o << ")";
245 }
246 }
247 else enregistrer_ccf(o,version);
248 o << ");" << std::endl;
249 }
250
251
252 void MG_ARETE::get_intersection(double t1, double t2,double* point_iners)
253 {
254 OT_VECTEUR_3D P1_INTERSECTION;
255
256 double U,W;
257
258 double xyz2[3];
259 double xyz1[3];
260 double drv1[3];
261 double drv2[3];
262
263 evaluer(t1,xyz1);
264 deriver(t1,drv1);
265
266 evaluer(t2,xyz2);
267 deriver(t2,drv2);
268
269 OT_VECTEUR_3D a(xyz1) ;
270 OT_VECTEUR_3D b(drv1) ;
271 OT_VECTEUR_3D c(xyz2) ;
272 OT_VECTEUR_3D d(drv2) ;
273 OT_VECTEUR_3D f=c-a;
274 d=-d;
275 double det1=b[0]*d[1]-b[1]*d[0];
276 double det2=b[0]*d[2]-b[2]*d[0];
277 double det3=b[1]*d[2]-b[2]*d[1];
278
279 bool c_fait=false;
280
281 if (fabs(det1)>=1e-10)
282 {
283 double u,w;
284 u=(d[1]*f[0]-d[0]*f[1])/det1;
285 w=(b[0]*f[1]-b[1]*f[0])/det1;
286 double verif_x=a[2] +u*b[2];
287 double verif_y=c[2]-w*d[2];
288 if (fabs(verif_x-verif_y)<0.000001)
289 {
290 U=u ;
291 W=w ;
292 c_fait=true;
293 }
294
295 }
296
297 if (!c_fait&&fabs(det2)>=1e-10)
298 {
299 double u,w;
300 u=(d[2]*f[0]-d[0]*f[2])/det2;
301 w=(b[0]*f[2]-b[2]*f[0])/det2;
302 double verif_x=a[1]+u*b[1];
303 double verif_y=c[1]-w*d[1];
304 if (fabs(verif_x-verif_y)<0.000001)
305 {
306 U=u ;
307 W=w ;
308 c_fait=true;
309 }
310 }
311
312 if (!c_fait&&fabs(det3)>=1e-10)
313 {
314 double u,w;
315 u=(d[2]*f[1]-d[1]*f[2])/det3;
316 w=(b[1]*f[2]-b[2]*f[1])/det3;
317 double verif_x=a[0] +u*b[0];
318 double verif_y=c[0]-w*d[0];
319 if (fabs(verif_x-verif_y)<0.000001)
320 {
321 U=u ;
322 W=w ;
323 c_fait=true;
324 }
325 }
326
327 P1_INTERSECTION=a+U*b;
328
329 point_iners[0]=P1_INTERSECTION[0];
330 point_iners[1]=P1_INTERSECTION[1];
331 point_iners[2]=P1_INTERSECTION[2];
332
333 }
334
335
336
337 void MG_ARETE:: get_param_NURBS(int& indx_premier_ptctr,TPL_LISTE_ENTITE<double> &param)
338 {
339
340 TPL_LISTE_ENTITE<double> param1;
341 int type=courbe->get_type_geometrique(param1);
342
343 double tcs1=get_tmin();
344 double tcs2=get_tmax();
345 double angle=tcs1-tcs2;
346
347 if ((type==MGCo_CIRCLE||type==MGCo_ELLIPSE)&&(tcs1!=tcs2)&&fabs(angle)!=2*PI)
348 {
349
350 double xyz[3];
351 param.ajouter(1);
352 param.ajouter(4);
353 param.ajouter(0);
354 param.ajouter(7);
355 param.ajouter(0);
356
357 double arc_lengh= get_longueur(tcs1,tcs2);
358 double rayon=arc_lengh/fabs(angle);
359
360 double dif=fabs(fabs(angle)-PI);
361 if (fabs(fabs(angle)-PI)<=1e-6||fabs(angle)-PI<0.) //ouverture inf ou egale a PI
362 {
363
364
365 double request_lengh=arc_lengh/2;
366
367 param.ajouter(0);
368 param.ajouter(0);
369 param.ajouter(0);
370 param.ajouter(0.5);
371 param.ajouter(0.5);
372 param.ajouter(1);
373 param.ajouter(1);
374 param.ajouter(1);
375
376
377 double xyz1[3];
378 double xyz1_inters[3];
379 double xyz2_inters[3];
380 double xyz2[3];
381 double xyz3[3];
382 double t_mid=get_t(tcs1,tcs2,request_lengh);
383
384 evaluer(tcs1,xyz1);
385
386 param.ajouter(xyz1[0]);
387 param.ajouter(xyz1[1]);
388 param.ajouter(xyz1[2]);
389 param.ajouter(1);
390
391 get_intersection(tcs1, t_mid,xyz1_inters);
392
393 param.ajouter(xyz1_inters[0]);
394 param.ajouter(xyz1_inters[1]);
395 param.ajouter(xyz1_inters[2]);
396 param.ajouter(0.5);
397
398 get_intersection(t_mid, tcs2,xyz2_inters);
399
400 param.ajouter(xyz2_inters[0]);
401 param.ajouter(xyz2_inters[1]);
402 param.ajouter(xyz2_inters[2]);
403 param.ajouter(0.5);
404
405 evaluer(tcs2,xyz3);
406
407 param.ajouter(xyz3[0]);
408 param.ajouter(xyz3[1]);
409 param.ajouter(xyz3[2]);
410 param.ajouter(1);
411
412 indx_premier_ptctr=13;
413
414 }
415
416 if (fabs(fabs(angle)-PI)>1e-6&&fabs(angle)-PI>0.)
417 {
418
419 double xyz[3];
420
421 double lengh_at_mid=arc_lengh;
422 double request_lengh=lengh_at_mid/2.;
423 double tmid=get_t(tcs1,tcs2,request_lengh);
424 double tmid1=get_t(tcs1,tmid,request_lengh/2.);
425 double tmid2=get_t(tmid,tcs2,request_lengh/2.);
426
427 param.ajouter(0);
428 param.ajouter(0);
429 param.ajouter(0);
430 param.ajouter(0.25);
431 param.ajouter(0.25);
432 param.ajouter(0.5);
433 param.ajouter(0.5);
434 param.ajouter(0.75);
435 param.ajouter(0.75);
436 param.ajouter(1);
437 param.ajouter(1);
438 param.ajouter(1);
439
440
441
442 evaluer(tcs1,xyz);
443
444 param.ajouter(xyz[0]);
445 param.ajouter(xyz[1]);
446 param.ajouter(xyz[2]);
447 param.ajouter(1);
448
449
450 get_intersection(tcs1,tmid1,xyz);
451
452 param.ajouter(xyz[0]);
453 param.ajouter(xyz[1]);
454 param.ajouter(xyz[2]);
455 param.ajouter(0.5);
456
457 get_intersection(tmid1, tmid,xyz);
458
459 param.ajouter(xyz[0]);
460 param.ajouter(xyz[1]);
461 param.ajouter(xyz[2]);
462 param.ajouter(0.5);
463
464 evaluer(tmid,xyz);
465
466 param.ajouter(xyz[0]);
467 param.ajouter(xyz[1]);
468 param.ajouter(xyz[2]);
469 param.ajouter(1);
470
471
472 get_intersection(tmid,tmid2,xyz);
473
474 param.ajouter(xyz[0]);
475 param.ajouter(xyz[1]);
476 param.ajouter(xyz[2]);
477 param.ajouter(0.5);
478
479
480 get_intersection(tmid2,tcs2,xyz);
481
482 param.ajouter(xyz[0]);
483 param.ajouter(xyz[1]);
484 param.ajouter(xyz[2]);
485 param.ajouter(0.5);
486
487 evaluer(tcs2,xyz);
488
489 param.ajouter(xyz[0]);
490 param.ajouter(xyz[1]);
491 param.ajouter(xyz[2]);
492 param.ajouter(1);
493
494 indx_premier_ptctr=17;
495 }
496
497 }
498
499
500 if ((type==MGCo_ELLIPSE||type==MGCo_CIRCLE)&&OPERATEUR::egal(fabs(angle),2*PI,0.000001))
501 courbe->get_param_NURBS(indx_premier_ptctr,param);
502
503 if (type==MGCo_BSPLINE)
504 courbe->get_param_NURBS(indx_premier_ptctr,param);
505
506 if (type==MGCo_LINE)
507 {
508
509 param.ajouter(1);
510 param.ajouter(2);
511 param.ajouter(0);
512 param.ajouter(2);
513 param.ajouter(0);
514
515
516 param.ajouter(0);
517 param.ajouter(0);
518 param.ajouter(1);
519 param.ajouter(1);
520
521 double xyz1[3];
522 double xyz2[3];
523
524 double tcs1=get_tmin();
525 double tcs2=get_tmax();
526
527 evaluer(tcs1,xyz1);
528
529 param.ajouter(xyz1[0]);
530 param.ajouter(xyz1[1]);
531 param.ajouter(xyz1[2]);
532 param.ajouter(1);
533
534 evaluer(tcs2,xyz2);
535
536 param.ajouter(xyz2[0]);
537 param.ajouter(xyz2[1]);
538 param.ajouter(xyz2[2]);
539 param.ajouter(1);
540
541 indx_premier_ptctr=9;
542 }
543
544 }
545
546 struct integrale {
547 double ti;
548 double li;
549 };
550
551 double MG_ARETE:: get_t(double t1,double t2,double longueur_voulue,double dt)
552 {
553 double pt[3];
554 double t,ti,tii;
555
556 //double longueur_voulue,longueur=0;
557 double longueur=0;
558 std::vector<integrale> tab_integrale;
559 //longueur_voulue= alpha *get_longueur(t1,t2);
560
561 integrale pas;
562 ti=t1;
563 do
564 {
565 tii=ti+(t2-t1)*dt;
566
567 t=0.7886751345*ti+0.2113248654*tii;
568 deriver(t,pt);
569
570 double facteur=pow(pt[0],2)+pow(pt[1],2)+pow(pt[2],2);
571 longueur=longueur+0.5*(tii-ti)*sqrt(facteur);
572 t=0.7886751345*tii+0.2113248654*ti;
573 deriver(t,pt);
574 facteur=pow(pt[0],2)+pow(pt[1],2)+pow(pt[2],2);
575 longueur=longueur+0.5*(tii-ti)*sqrt(facteur);
576 integrale pas;
577 pas.ti=ti;
578 pas.li=longueur;
579 tab_integrale.insert(tab_integrale.end(),pas);
580 ti=tii;
581 }
582 while (fabs(tab_integrale.back().li)<fabs(longueur_voulue));
583
584 //double t_trouve=0.7886751345*tab_integrale.back().ti+0.2113248654*tab_integrale.back().ti+(t2-t1)*dt;
585 double tn_1=tab_integrale.back().ti;
586 double ln_1=tab_integrale[tab_integrale.size()-2].li;
587 double tn=ti;
588 double ln=tab_integrale.back().li;
589
590
591 double t_trouve=tn_1+(tn-tn_1)/(ln-ln_1)*(longueur_voulue-ln_1);
592 return t_trouve ;
593 }
594
595 void MG_ARETE::recupere_resultat(class FEM_SOLUTION* sol,int numchamps,char* fichier)
596 {
597 FEM_MAILLAGE *mai=sol->get_maillage();
598 MG_SOMMET* som1=cosommet1->get_sommet();
599 MG_SOMMET* som2=cosommet2->get_sommet();
600 sol->active_solution(numchamps);
601
602 FEM_NOEUD* noeudcourant=NULL;
603 for (int i=0;i<som1->get_lien_fem_maillage()->get_nb();i++)
604 {
605 FEM_NOEUD* no=(FEM_NOEUD*)som1->get_lien_fem_maillage()->get(i);
606 if (mai->get_fem_noeudid(no->get_id())==no) noeudcourant=no;
607 }
608 FEM_NOEUD* noeudfin=NULL;
609 for (int i=0;i<som2->get_lien_fem_maillage()->get_nb();i++)
610 {
611 FEM_NOEUD* no=(FEM_NOEUD*)som2->get_lien_fem_maillage()->get(i);
612 if (mai->get_fem_noeudid(no->get_id())==no) noeudfin=no;
613 }
614 int n1=0;
615 int n2;
616 if (mai->get_degre()==1) n2=1;
617 if (mai->get_degre()==2) n2=2;
618 double *xyz=noeudcourant->get_coord();
619 double t;
620 inverser(t,xyz);
621 std::vector<double> resultat;
622 double tprec=t;
623 resultat.push_back(t);
624 resultat.push_back(0.);
625 resultat.push_back(noeudcourant->get_solution());
626 while (noeudcourant!=noeudfin)
627 {
628 FEM_ELEMENT1 *segcourant=NULL;
629 for (int i=0;i<get_lien_fem_maillage()->get_nb();i++)
630 {
631 FEM_ELEMENT1* seg=(FEM_ELEMENT1*)get_lien_fem_maillage()->get(i);
632 if (mai->get_fem_element1id(seg->get_id())==seg)
633 {
634 if ((seg->get_fem_noeud(n1)==noeudcourant)||(seg->get_fem_noeud(n2)==noeudcourant))
635 if (seg!=segcourant)
636 segcourant=seg;
637 }
638 }
639 FEM_NOEUD* autrenoeud=NULL;
640 if (segcourant->get_fem_noeud(n1)==noeudcourant) autrenoeud=segcourant->get_fem_noeud(n2);
641 if (segcourant->get_fem_noeud(n2)==noeudcourant) autrenoeud=segcourant->get_fem_noeud(n1);
642 FEM_NOEUD *noprec=noeudcourant;
643 for (int i=1;i<segcourant->get_nb_fem_noeud();i++)
644 {
645 FEM_NOEUD* no=segcourant->get_fem_noeud(i);
646 double *xyz=no->get_coord();
647 double t2;
648 inverser(t2,xyz);
649 if (get_courbe()->est_periodique())
650 if (t2<tprec) t2=t2+get_courbe()->get_periode();
651 double s=resultat[resultat.size()-2];
652 int pas=32;
653 for (int j=0;j<pas;j++)
654 {
655 double ti=tprec+1.0*j*(t2-tprec)/pas;
656 double tii=tprec+1.0*(j+1)*(t2-tprec)/pas;
657 double tgauss1=0.7886751345*tii+0.2113248654*ti;
658 double dxyz[3];
659 deriver(tgauss1,dxyz);
660 s=s+0.5*(tii-ti)*sqrt(dxyz[0]*dxyz[0]+dxyz[1]*dxyz[1]+dxyz[2]*dxyz[2]);
661 double tgauss2=0.7886751345*ti+0.2113248654*tii;
662 deriver(tgauss2,dxyz);
663 s=s+0.5*(tii-ti)*sqrt(dxyz[0]*dxyz[0]+dxyz[1]*dxyz[1]+dxyz[2]*dxyz[2]);
664 }
665 resultat.push_back(t2);
666 resultat.push_back(s);
667 resultat.push_back(no->get_solution());
668 noprec=no;
669 tprec=t2;
670 }
671 noeudcourant=autrenoeud;
672 }
673 FILE *in=fopen(fichier,"wt");
674 char mess[500];
675 sprintf(mess,"%s_%s",(char*)sol->get_nom().c_str(),(char*)sol->get_legende(numchamps).c_str());
676 fprintf(in,"t;s;%s;\n",mess);
677 for (int i=0;i<resultat.size();i=i+3)
678 fprintf(in,"%lf;%lf;%lf;\n",resultat[i],resultat[i+1],resultat[i+2]);
679 fclose(in);
680 }