ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mtu/src/mg_arete.cpp
Revision: 1158
Committed: Thu Jun 13 22:18:49 2024 UTC (11 months, 1 week ago) by francois
File size: 18330 byte(s)
Log Message:
compatibilité Ubuntu 22.04
Suppression des refeences à Windows
Ajout d'une banière

File Contents

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