ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/mg_arete.cpp
Revision: 763
Committed: Wed Dec 2 19:55:53 2015 UTC (9 years, 5 months ago) by francois
File size: 17563 byte(s)
Log Message: 
Le fichier MAGiC est maintenant versionné. LA version actuelle est 2.0. L'ancienne version est 1.0.
Tout est transparent pour l'utilisateur. Les vieilles versions sont lisibles mais les nouveaux enregistrements sont dans la version la plus récente.
Changement des conditions aux limites : ajout d'un parametre pour dire si la condition numerique est une valeur ou une formule ou un lien vers une autre entité magic.
Les parametres pour saisir sont maintenant -ccf -ccfi -ccff -ccft -ccfit -ccfft

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