ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/outil/src/ot_mathematique.cpp
Revision: 86
Committed: Fri May 2 14:26:38 2008 UTC (17 years ago) by souaissa
Original Path: magic/lib/outil/outil/src/ot_mathematique.cpp
File size: 25206 byte(s)
Log Message: 
Mise à jour des classes vct;
Mise à jour de la classe ot_tenseur;

File Contents

# User Rev Content
1 souaissa 58 //------------------------------------------------------------
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     // ot_mathematique.cpp
16     //
17     //------------------------------------------------------------
18     //------------------------------------------------------------
19     // COPYRIGHT 2000
20     // Version du 02/03/2006 à 11H23
21     //------------------------------------------------------------
22     //------------------------------------------------------------
23    
24    
25     #include "gestionversion.h"
26 souaissa 70 #include <iomanip.h>
27 souaissa 58 #include <math.h>
28 souaissa 86 #include <assert.h>
29 souaissa 58 #include <algorithm>
30     #include "ot_mathematique.h"
31    
32     using namespace std;
33    
34    
35     OT_VECTEUR_3D::OT_VECTEUR_3D(double x,double y,double z)
36     {
37     valeur[0]=x;
38     valeur[1]=y;
39     valeur[2]=z;
40     }
41    
42    
43     OT_VECTEUR_3D::OT_VECTEUR_3D(double *xyz)
44     {
45     valeur[0]=xyz[0];
46     valeur[1]=xyz[1];
47     valeur[2]=xyz[2];
48     }
49    
50     OT_VECTEUR_3D::OT_VECTEUR_3D(double *xyz1,double *xyz2)
51     {
52     valeur[0]=xyz2[0]-xyz1[0];
53     valeur[1]=xyz2[1]-xyz1[1];
54     valeur[2]=xyz2[2]-xyz1[2];
55     }
56    
57     OT_VECTEUR_3D::OT_VECTEUR_3D(const OT_VECTEUR_3D& mdd)
58     {
59     valeur[0]=mdd.valeur[0];
60     valeur[1]=mdd.valeur[1];
61     valeur[2]=mdd.valeur[2];
62     }
63    
64     OT_VECTEUR_3D::OT_VECTEUR_3D(OT_VECTEUR_3D& mdd)
65     {
66     valeur[0]=mdd.valeur[0];
67     valeur[1]=mdd.valeur[1];
68     valeur[2]=mdd.valeur[2];
69     }
70    
71     OT_VECTEUR_3D::OT_VECTEUR_3D(void)
72     {
73     valeur[0]=0.;
74     valeur[1]=0.;
75     valeur[2]=0.;
76     }
77    
78     double OT_VECTEUR_3D::get_x(void) const
79     {
80     return valeur[0];
81     }
82    
83    
84     double OT_VECTEUR_3D::get_y(void) const
85     {
86     return valeur[1];
87     }
88    
89    
90     double OT_VECTEUR_3D::get_z(void) const
91     {
92     return valeur[2];
93     }
94    
95     void OT_VECTEUR_3D::change_x(double x)
96     {
97     valeur[0]=x;
98     }
99    
100     void OT_VECTEUR_3D::change_y(double y)
101     {
102     valeur[1]=y;
103     }
104    
105     void OT_VECTEUR_3D::change_z(double z)
106     {
107     valeur[2]=z;
108     }
109    
110     double* OT_VECTEUR_3D::get_xyz(void)
111     {
112     return valeur;
113     }
114    
115     OT_VECTEUR_3D::operator const double* () const
116     {
117     return valeur;
118     }
119    
120     OT_VECTEUR_3D::operator double* ()
121     {
122     return valeur;
123     }
124    
125     double OT_VECTEUR_3D::operator[] (int i) const
126     {
127     return valeur[i];
128     }
129    
130     double & OT_VECTEUR_3D::operator[] (int i)
131     {
132     return valeur[i];
133     }
134    
135     double OT_VECTEUR_3D::operator() (int i) const
136     {
137     return valeur[i];
138     }
139     double & OT_VECTEUR_3D::operator() (int i)
140     {
141     return valeur[i];
142     }
143    
144     double OT_VECTEUR_3D::get_longueur(void) const
145     {
146     return sqrt(valeur[0]*valeur[0]+valeur[1]*valeur[1]+valeur[2]*valeur[2]);
147     }
148    
149     double OT_VECTEUR_3D::get_longueur2(void) const
150     {
151     return (valeur[0]*valeur[0]+valeur[1]*valeur[1]+valeur[2]*valeur[2]);
152     }
153    
154    
155     void OT_VECTEUR_3D::norme(void)
156     {
157     double longueur=get_longueur();
158     if (longueur != 0)
159     {
160     valeur[0]=valeur[0]/longueur;
161     valeur[2]=valeur[2]/longueur;
162     valeur[1]=valeur[1]/longueur;
163     }
164     }
165    
166     OT_VECTEUR_3D OT_VECTEUR_3D::gram_shmidt(const OT_VECTEUR_3D& vint)
167     {
168     OT_VECTEUR_3D v(0.,0.,0.);
169     double alpha=(*this)*vint;
170     v= *this-alpha*vint;
171     return v;
172     }
173    
174     OT_VECTEUR_3D OT_VECTEUR_3D::unite(int i)
175     {
176     OT_VECTEUR_3D v(0.,0.,0.);
177     v[i]=1.;
178     return v;
179     }
180    
181    
182     int OT_VECTEUR_3D::compare_valeur (const OT_VECTEUR_3D& mdd) const
183     {
184     return memcmp(valeur, mdd.valeur, 3*sizeof(double));
185     }
186    
187     bool OT_VECTEUR_3D::operator== (const OT_VECTEUR_3D& mdd) const
188     {
189     return compare_valeur(mdd) == 0;
190     }
191    
192     bool OT_VECTEUR_3D::operator!= (const OT_VECTEUR_3D& mdd) const
193     {
194     return compare_valeur(mdd) != 0;
195     }
196    
197     bool OT_VECTEUR_3D::operator< (const OT_VECTEUR_3D& mdd) const
198     {
199     return compare_valeur(mdd) < 0;
200     }
201    
202     bool OT_VECTEUR_3D::operator<= (const OT_VECTEUR_3D& mdd) const
203     {
204     return compare_valeur(mdd) <= 0;
205     }
206    
207     bool OT_VECTEUR_3D::operator> (const OT_VECTEUR_3D& mdd) const
208     {
209     return compare_valeur(mdd) > 0;
210     }
211    
212     bool OT_VECTEUR_3D::operator>= (const OT_VECTEUR_3D& mdd) const
213     {
214     return compare_valeur(mdd) >= 0;
215     }
216    
217     std::ostream & operator << ( std::ostream & __os, const OT_VECTEUR_3D & __vec)
218     {
219     __os << "{ "<< __vec[0]<<", "<<__vec[1]<<", "<<__vec[2]<<"}, ";
220     return __os;
221     }
222    
223     OT_VECTEUR_3D operator/ (const OT_VECTEUR_3D& mdd1, double diviseur)
224     {
225     OT_VECTEUR_3D kQuot;
226    
227     if ( diviseur != (double)0.0 )
228     {
229     double facteur = ((double)1.0)/diviseur;
230     kQuot.valeur[0] = facteur*mdd1.valeur[0];
231     kQuot.valeur[1] = facteur*mdd1.valeur[1];
232     kQuot.valeur[2] = facteur*mdd1.valeur[2];
233     }
234     else
235     {
236     kQuot.valeur[0] = 1E300;
237     kQuot.valeur[1] = 1E300;
238     kQuot.valeur[2] = 1E300;
239     }
240    
241     return kQuot;
242     }
243    
244     OT_VECTEUR_3D operator- (const OT_VECTEUR_3D& mdd1)
245     {
246     return OT_VECTEUR_3D(
247     -mdd1.valeur[0],
248     -mdd1.valeur[1],
249     -mdd1.valeur[2]);
250     }
251    
252     double operator*(const OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
253     {
254     double tmp;
255     tmp=mdd1.valeur[0]*mdd2.valeur[0]+mdd1.valeur[1]*mdd2.valeur[1]+mdd1.valeur[2]*mdd2.valeur[2];
256     return tmp;
257     }
258    
259     OT_VECTEUR_3D operator*(double mdd1,const OT_VECTEUR_3D& mdd2)
260     {
261     OT_VECTEUR_3D tmp;
262     tmp.valeur[0]=mdd1*mdd2.valeur[0];
263     tmp.valeur[1]=mdd1*mdd2.valeur[1];
264     tmp.valeur[2]=mdd1*mdd2.valeur[2];
265     return tmp;
266     }
267    
268     OT_VECTEUR_3D& OT_VECTEUR_3D::operator+= (const OT_VECTEUR_3D& mdd)
269     {
270     valeur[0] += mdd.valeur[0];
271     valeur[1] += mdd.valeur[1];
272     valeur[2] += mdd.valeur[2];
273     return *this;
274     }
275     OT_VECTEUR_3D& OT_VECTEUR_3D::operator-= (const OT_VECTEUR_3D& mdd)
276     {
277     valeur[0] -= mdd.valeur[0];
278     valeur[1] -= mdd.valeur[1];
279     valeur[2] -= mdd.valeur[2];
280     return *this;
281     }
282     OT_VECTEUR_3D& OT_VECTEUR_3D::operator*= (double facteur)
283     {
284     valeur[0] *= facteur;
285     valeur[1] *= facteur;
286     valeur[2] *= facteur;
287     return *this;
288     }
289     OT_VECTEUR_3D& OT_VECTEUR_3D::operator/= (double diviseur)
290     {
291     if ( diviseur != (double)0.0 )
292     {
293     double facteur = ((double)1.0)/diviseur;
294     valeur[0] *= facteur;
295     valeur[1] *= facteur;
296     valeur[2] *= facteur;
297     }
298     else
299     {
300     valeur[0] = 1E300;
301     valeur[1] = 1E300;
302     valeur[2] = 1E300;
303     }
304    
305     return *this;
306     }
307    
308     OT_VECTEUR_3D operator&(const OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
309     {
310     OT_VECTEUR_3D tmp;
311     tmp.valeur[0]=mdd1.valeur[1]*mdd2.valeur[2]-mdd1.valeur[2]*mdd2.valeur[1];
312     tmp.valeur[1]=mdd1.valeur[2]*mdd2.valeur[0]-mdd1.valeur[0]*mdd2.valeur[2];
313     tmp.valeur[2]=mdd1.valeur[0]*mdd2.valeur[1]-mdd1.valeur[1]*mdd2.valeur[0];
314     return tmp;
315     }
316    
317     double OT_VECTEUR_3D::diff(void)
318     {
319     return 0.33333333333333*(fabs(valeur[0])+fabs(valeur[1])+fabs(valeur[2]));
320     }
321    
322    
323     OT_MATRICE_3D::OT_MATRICE_3D(OT_VECTEUR_3D& colonne1,OT_VECTEUR_3D& colonne2,OT_VECTEUR_3D& colonne3)
324     {
325     vec[0]=colonne1;
326     vec[1]=colonne2;
327     vec[2]=colonne3;
328     }
329    
330     OT_MATRICE_3D::OT_MATRICE_3D(const OT_MATRICE_3D& mdd)
331     {
332     vec[0]=mdd.vec[0];
333     vec[1]=mdd.vec[1];
334     vec[2]=mdd.vec[2];
335     }
336    
337    
338     std::ostream & operator << ( std::ostream & __os, const OT_MATRICE_3D & __mat)
339     {
340    
341     __os << "{ ";
342     for (unsigned i=0; i<3; i++)
343     {
344     __os << "[ "<< __mat(i,0)<<" ]\t [ "<<__mat(i,1) <<" ]\t [ "<<__mat(i,2);
345     if (i+1 < 3)
346     __os << std::endl << " ";
347     }
348     __os << " }" << std::endl;
349     return __os;
350     }
351    
352     OT_VECTEUR_3D & OT_MATRICE_3D::operator [](int i)
353     {
354     return vec[i];
355     }
356    
357     OT_MATRICE_3D::OT_MATRICE_3D(double* t)
358     {
359     for(int i=0;i<3;i++)
360     {
361     for(int j=0;j<3;j++)
362     {
363     vec[i][j]=t[i*3+j];
364     }
365    
366     }
367     }
368    
369     double OT_MATRICE_3D::valeur(int iLigne, int iCol) const
370     {
371     return vec[iCol][iLigne];
372     }
373     double & OT_MATRICE_3D::valeur(int iLigne, int iCol)
374     {
375     return vec[iCol][iLigne];
376     }
377     double OT_MATRICE_3D::operator() (int iLigne, int iCol) const
378     {
379     return vec[iCol][iLigne];
380     }
381     double & OT_MATRICE_3D::operator() (int iLigne, int iCol)
382     {
383     return vec[iCol][iLigne];
384     }
385     double OT_MATRICE_3D::get_determinant()
386     {
387     double coFacteur0 = valeur(1,1)*valeur(2,2) - valeur(1,2)*valeur(2,1);
388     double coFacteur1 = valeur(1,2)*valeur(2,0) - valeur(1,0)*valeur(2,2);
389     double coFacteur2 = valeur(1,0)*valeur(2,1) - valeur(1,1)*valeur(2,0);
390     double determinant = valeur(0,0)*coFacteur0 + valeur(0,1)*coFacteur1 + valeur(0,2)*coFacteur2;
391     return determinant;
392     }
393     void OT_MATRICE_3D::identite(void)
394     {
395     for (unsigned i = 0; i<3; i++)
396     for (unsigned j = 0; j < 3; j++)
397     vec[i][j] = (i==j) ? 1.0 : 0.0 ;
398     }
399    
400     void OT_MATRICE_3D::transpose(OT_MATRICE_3D& res) const
401     {
402     for (unsigned i = 0; i<3; i++)
403     for (unsigned j = 0; j < 3; j++)
404     res.vec[i][j] = vec[j][i];
405     }
406    
407     OT_MATRICE_3D OT_MATRICE_3D::inverse () const
408     {
409     // Invert a 3x3 using cofactors. This is faster than using a generic
410     // Gaussian elimination because of the loop overhead of such a method.
411    
412     OT_MATRICE_3D matInverse;
413    
414     matInverse.valeur(0,0) = valeur(1,1)*valeur(2,2) - valeur(1,2)*valeur(2,1);
415     matInverse.valeur(0,1) = valeur(0,2)*valeur(2,1) - valeur(0,1)*valeur(2,2);
416     matInverse.valeur(0,2) = valeur(0,1)*valeur(1,2) - valeur(0,2)*valeur(1,1);
417     matInverse.valeur(1,0) = valeur(1,2)*valeur(2,0) - valeur(1,0)*valeur(2,2);
418     matInverse.valeur(1,1) = valeur(0,0)*valeur(2,2) - valeur(0,2)*valeur(2,0);
419     matInverse.valeur(1,2) = valeur(0,2)*valeur(1,0) - valeur(0,0)*valeur(1,2);
420     matInverse.valeur(2,0) = valeur(1,0)*valeur(2,1) - valeur(1,1)*valeur(2,0);
421     matInverse.valeur(2,1) = valeur(0,1)*valeur(2,0) - valeur(0,0)*valeur(2,1);
422     matInverse.valeur(2,2) = valeur(0,0)*valeur(1,1) - valeur(0,1)*valeur(1,0);
423    
424     double det = valeur(0,0)*matInverse(0,0) + valeur(0,1)*matInverse(1,0)+
425     valeur(0,2)*matInverse(2,0);
426    
427     if ( fabs(det) <= 1E-100 )
428     {
429     matInverse.identite();
430     return matInverse;
431     }
432    
433     double invDet = 1/det;
434    
435     for (int i=0; i<3; i++)
436     for (int j=0; j<3; j++)
437     matInverse.vec[i][j] *= invDet;
438    
439     return matInverse;
440     }
441    
442     void OT_MATRICE_3D::change_vecteur1(OT_VECTEUR_3D v)
443     {
444     vec[0]=v;
445     }
446     void OT_MATRICE_3D::change_vecteur2(OT_VECTEUR_3D v)
447     {
448     vec[1]=v;
449     }
450     void OT_MATRICE_3D::change_vecteur3(OT_VECTEUR_3D v)
451     {
452     vec[2]=v;
453     }
454     OT_VECTEUR_3D& OT_MATRICE_3D::get_vecteur1(void)
455     {
456     return vec[0];
457     }
458     OT_VECTEUR_3D& OT_MATRICE_3D::get_vecteur2(void)
459     {
460     return vec[1];
461     }
462     OT_VECTEUR_3D& OT_MATRICE_3D::get_vecteur3(void)
463     {
464     return vec[2];
465     }
466    
467    
468    
469 souaissa 69 OT_MATRICE_3D operator+(const OT_MATRICE_3D& mdd1,const OT_MATRICE_3D& mdd2)
470 foucault 27 {
471 souaissa 58 OT_MATRICE_3D tmp;
472     for(int i=0;i<3;i++) {
473     for(int j=0;j<3;j++) tmp(i,j)= mdd1(i,j)+mdd2(i,j);
474     }
475     return tmp;
476 foucault 27 }
477 souaissa 58
478    
479     OT_VECTEUR_3D operator*(const OT_MATRICE_3D& mdd1,const OT_VECTEUR_3D& mdd2)
480     {
481     OT_VECTEUR_3D tmp;
482     tmp(0)=mdd1(0,0)*mdd2(0)+mdd1(0,1)*mdd2(1)+mdd1(0,2)*mdd2(2);
483     tmp(1)=mdd1(1,0)*mdd2(0)+mdd1(1,1)*mdd2(1)+mdd1(1,2)*mdd2(2);
484     tmp(2)=mdd1(2,0)*mdd2(0)+mdd1(2,1)*mdd2(1)+mdd1(2,2)*mdd2(2);
485     return tmp;}
486    
487     OT_VECTEUR_3D operator+(const OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
488     {
489     OT_VECTEUR_3D tmp;
490     tmp(0)=mdd1(0)+mdd2(0);
491     tmp(1)=mdd1(1)+mdd2(1);
492     tmp(2)=mdd1(2)+mdd2(2);
493     return tmp;
494     }
495    
496     OT_VECTEUR_3D operator-(const class OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
497     {
498     OT_VECTEUR_3D tmp;
499     tmp.valeur[0]=mdd1.valeur[0]-mdd2.valeur[0];
500     tmp.valeur[1]=mdd1.valeur[1]-mdd2.valeur[1];
501     tmp.valeur[2]=mdd1.valeur[2]-mdd2.valeur[2];
502     return tmp;
503     }
504    
505     OT_MATRICE_3D operator*(const OT_MATRICE_3D& mdd1,const OT_MATRICE_3D& mdd2)
506     {
507     OT_MATRICE_3D tmp;
508     for (int i=0; i<3; i++)
509     for (int j=0; j<3; j++)
510     {
511     tmp(i,j) = 0;
512     for (int k=0; k<3; k++)
513     tmp(i,j) += mdd1(i,k)*mdd2(k,j);
514     }
515     return tmp;
516     }
517    
518     #ifdef ot_vecteur_4d
519     OT_VECTEUR_4D::OT_VECTEUR_4D (double v[4])
520     {
521     for (int i=0; i<4; i++)
522     valeur[i]=v[i];
523     }
524    
525     OT_VECTEUR_4D::OT_VECTEUR_4D (double x, double y, double z, double w)
526     {
527     valeur[0] = x;
528     valeur[1] = y;
529     valeur[2] = z;
530     valeur[3] = w;
531     }
532    
533     OT_VECTEUR_4D::OT_VECTEUR_4D (const OT_VECTEUR_4D & v)
534     {
535     for (int i=0; i<4; i++)
536     valeur[i]=v[i];
537     }
538    
539    
540     double & OT_VECTEUR_4D::operator()(int i)
541     {
542     return valeur[i];
543     }
544    
545     double OT_VECTEUR_4D::operator()(int i) const
546     {
547     return valeur[i];
548     }
549    
550     double & OT_VECTEUR_4D::operator[](int i)
551     {
552     return valeur[i];
553     }
554    
555     double OT_VECTEUR_4D::operator[](int i) const
556     {
557     return valeur[i];
558     }
559    
560     OT_VECTEUR_4D::operator const double* () const
561     {
562     return valeur;
563     }
564    
565     OT_VECTEUR_4D::operator double* ()
566     {
567     return valeur;
568     }
569    
570     double OT_VECTEUR_4D::get_x(void) const
571     { return valeur[0]; }
572     double & OT_VECTEUR_4D::x(void)
573     { return valeur[0]; }
574     double OT_VECTEUR_4D::get_y(void) const
575     { return valeur[1]; }
576     double & OT_VECTEUR_4D::y(void)
577     { return valeur[1]; }
578     double OT_VECTEUR_4D::get_z(void) const
579     { return valeur[2]; }
580     double & OT_VECTEUR_4D::z(void)
581     { return valeur[2]; }
582     double OT_VECTEUR_4D::get_w(void) const
583     { return valeur[3]; }
584     double & OT_VECTEUR_4D::w(void)
585     { return valeur[3]; }
586     void OT_VECTEUR_4D::change_x(double x)
587     { valeur[0] = x; }
588     void OT_VECTEUR_4D::change_y(double y)
589     { valeur[1] = y; }
590     void OT_VECTEUR_4D::change_z(double z)
591     { valeur[2] = z; }
592     void OT_VECTEUR_4D::change_w(double w)
593     { valeur[3] = w; }
594     double* OT_VECTEUR_4D::get_xyzw(void)
595     { return valeur; }
596    
597     OT_VECTEUR_4D OT_VECTEUR_4D::operator+ (const OT_VECTEUR_4D & a)
598     {
599     return OT_VECTEUR_4D(
600     a(0)+valeur[0],a(1)+valeur[1],a(2)+valeur[2],a(3)+valeur[3]
601     );
602     }
603    
604     OT_VECTEUR_4D OT_VECTEUR_4D::operator- (const OT_VECTEUR_4D & a)
605     {
606     return OT_VECTEUR_4D(
607     valeur[0]-a(0),valeur[1]-a(1),valeur[2]-a(2),valeur[3]-a(3)
608     );
609     }
610    
611     OT_VECTEUR_4D operator* (const OT_VECTEUR_4D & rkV, const double a)
612     {
613     return OT_VECTEUR_4D(
614     a*rkV[0],a*rkV[1],a*rkV[2],a*rkV[3]
615     );
616     }
617     OT_VECTEUR_4D operator* (const double a, const OT_VECTEUR_4D & rkV)
618     {
619     return OT_VECTEUR_4D(
620     a*rkV[0],a*rkV[1],a*rkV[2],a*rkV[3]
621     );
622     }
623    
624     OT_VECTEUR_4D operator/ (const OT_VECTEUR_4D & rkV, const double a)
625     {
626     double inv = 1/a;
627     return OT_VECTEUR_4D(
628     inv*rkV[0],inv*rkV[1],inv*rkV[2],inv*rkV[3]
629     );
630     }
631    
632    
633     OT_VECTEUR_4D operator/ (const double a, const OT_VECTEUR_4D & rkV)
634     {
635     double inv = 1/a;
636     return OT_VECTEUR_4D(
637     inv*rkV[0],inv*rkV[1],inv*rkV[2],inv*rkV[3]
638     );
639     }
640    
641     double OT_VECTEUR_4D::operator* (const OT_VECTEUR_4D & a)
642     {
643     return a[0]*valeur[0]+a[1]*valeur[1]+a[2]*valeur[2]+a[3]*valeur[3];
644     }
645    
646    
647     OT_VECTEUR_4D& OT_VECTEUR_4D::operator+= (const OT_VECTEUR_4D& rkV)
648     {
649     valeur[0] += rkV.valeur[0];
650     valeur[1] += rkV.valeur[1];
651     valeur[2] += rkV.valeur[2];
652     valeur[3] += rkV.valeur[3];
653     return *this;
654     }
655    
656     OT_VECTEUR_4D& OT_VECTEUR_4D::operator-= (const OT_VECTEUR_4D& rkV)
657     {
658     valeur[0] -= rkV.valeur[0];
659     valeur[1] -= rkV.valeur[1];
660     valeur[2] -= rkV.valeur[2];
661     valeur[3] -= rkV.valeur[3];
662     return *this;
663     }
664     OT_VECTEUR_4D& OT_VECTEUR_4D::operator*= (double fScalar)
665     {
666     valeur[0] *= fScalar;
667     valeur[1] *= fScalar;
668     valeur[2] *= fScalar;
669     valeur[3] *= fScalar;
670     return *this;
671     }
672     OT_VECTEUR_4D& OT_VECTEUR_4D::operator/= (double fScalar)
673     {
674     double a = 1/fScalar;
675     valeur[0] *= a;
676     valeur[1] *= a;
677     valeur[2] *= a;
678     valeur[3] *= a;
679     return *this;
680     }
681     OT_VECTEUR_4D &OT_VECTEUR_4D::operator= (double a)
682     {
683     valeur[0]=valeur[1]=valeur[2]=valeur[3]=a;
684     return *this;
685     }
686    
687     std::ostream & operator << ( std::ostream & __os, const OT_VECTEUR_4D & __vec)
688     {
689     __os << "{ "<< __vec[0]<<", "<<__vec[1]<<", "<<__vec[2]<<", "<<__vec[3]<<"}, ";
690     return __os;
691     }
692     #endif // ot_vecteur_4d
693    
694 souaissa 69
695    
696    
697     OT_VECTEUR_4DD::OT_VECTEUR_4DD (double2 v[4])
698     {
699     for (int i=0; i<4; i++)
700     valeur[i]=v[i];
701     }
702    
703     OT_VECTEUR_4DD::OT_VECTEUR_4DD (double2 x, double2 y, double2 z, double2 w)
704     {
705     valeur[0] = x;
706     valeur[1] = y;
707     valeur[2] = z;
708     valeur[3] = w;
709     }
710    
711     OT_VECTEUR_4DD::OT_VECTEUR_4DD (const OT_VECTEUR_4DD & v)
712     {
713     for (int i=0; i<4; i++)
714     valeur[i]=v[i];
715     }
716    
717    
718     double2 & OT_VECTEUR_4DD::operator()(int i)
719     {
720     return valeur[i];
721     }
722    
723     double2 OT_VECTEUR_4DD::operator()(int i) const
724     {
725     return valeur[i];
726     }
727    
728     double2 & OT_VECTEUR_4DD::operator[](int i)
729     {
730     return valeur[i];
731     }
732    
733     double2 OT_VECTEUR_4DD::operator[](int i) const
734     {
735     return valeur[i];
736     }
737    
738     OT_VECTEUR_4DD::operator const double2* () const
739     {
740     return valeur;
741     }
742    
743     OT_VECTEUR_4DD::operator double2* ()
744     {
745     return valeur;
746     }
747    
748     OT_VECTEUR_4DD& OT_VECTEUR_4DD::operator =(const OT_VECTEUR_4DD& v4d)
749     {
750     if (this==&v4d) return *this;
751     valeur [0]=v4d.get_x();
752     valeur [1]=v4d.get_y();
753     valeur [2]=v4d.get_z();
754     valeur [3]=v4d.get_w();
755     return *this;
756     }
757 souaissa 86 double2 OT_VECTEUR_4DD::norme()
758     {
759     double2 nor= valeur [0]*valeur [0]+valeur [1]*valeur [1]+valeur [2]*valeur [2]+valeur [3]*valeur [3];
760     nor=nor^0.5;
761     return nor;
762     }
763     OT_VECTEUR_4DD OT_VECTEUR_4DD::vecteur_norme()
764     {
765     OT_VECTEUR_4DD nor;
766     double2 zero=0.;
767     double2 norme= this->norme() ;
768     assert(norme!=zero) ;
769     nor=(*this)*(1./norme);
770     return nor;
771     }
772 souaissa 69
773     double2 OT_VECTEUR_4DD::get_x(void) const
774     { return valeur[0]; }
775     double2 & OT_VECTEUR_4DD::x(void)
776     { return valeur[0]; }
777     double2 OT_VECTEUR_4DD::get_y(void) const
778     { return valeur[1]; }
779     double2 & OT_VECTEUR_4DD::y(void)
780     { return valeur[1]; }
781     double2 OT_VECTEUR_4DD::get_z(void) const
782     { return valeur[2]; }
783     double2 & OT_VECTEUR_4DD::z(void)
784     { return valeur[2]; }
785     double2 OT_VECTEUR_4DD::get_w(void) const
786     { return valeur[3]; }
787     double2 & OT_VECTEUR_4DD::w(void)
788     { return valeur[3]; }
789     void OT_VECTEUR_4DD::change_x(double2 x)
790     { valeur[0] = x; }
791     void OT_VECTEUR_4DD::change_y(double2 y)
792     { valeur[1] = y; }
793     void OT_VECTEUR_4DD::change_z(double2 z)
794     { valeur[2] = z; }
795     void OT_VECTEUR_4DD::change_w(double2 w)
796     { valeur[3] = w; }
797     double2* OT_VECTEUR_4DD::get_xyzw(void)
798     { return valeur; }
799    
800     OT_VECTEUR_4DD OT_VECTEUR_4DD::operator+ (const OT_VECTEUR_4DD & a)
801     {
802     return OT_VECTEUR_4DD(
803     a(0)+valeur[0],a(1)+valeur[1],a(2)+valeur[2],a(3)+valeur[3]
804     );
805     }
806    
807     OT_VECTEUR_4DD OT_VECTEUR_4DD::operator- (const OT_VECTEUR_4DD & a)
808     {
809     return OT_VECTEUR_4DD(
810     valeur[0]-a(0),valeur[1]-a(1),valeur[2]-a(2),valeur[3]-a(3)
811     );
812     }
813    
814     OT_VECTEUR_4DD operator* (const OT_VECTEUR_4DD & rkV, const double2 a)
815     {
816     return OT_VECTEUR_4DD(
817     a*rkV[0],a*rkV[1],a*rkV[2],a*rkV[3]
818     );
819     }
820     OT_VECTEUR_4DD operator* (const double2 a, const OT_VECTEUR_4DD & rkV)
821     {
822     return OT_VECTEUR_4DD(
823     a*rkV[0],a*rkV[1],a*rkV[2],a*rkV[3]
824     );
825     }
826    
827     OT_VECTEUR_4DD operator/ (const OT_VECTEUR_4DD & rkV, const double2 a)
828     {
829     double2 inv = 1/a;
830     return OT_VECTEUR_4DD(
831     inv*rkV[0],inv*rkV[1],inv*rkV[2],inv*rkV[3]
832     );
833     }
834    
835    
836     OT_VECTEUR_4DD operator/ (const double2 a, const OT_VECTEUR_4DD & rkV)
837     {
838     double2 inv = 1/a;
839     return OT_VECTEUR_4DD(
840     inv*rkV[0],inv*rkV[1],inv*rkV[2],inv*rkV[3]
841     );
842     }
843 souaissa 86 OT_VECTEUR_4DD OT_VECTEUR_4DD::operator^(const OT_VECTEUR_4DD& v2)
844     {
845     OT_VECTEUR_4DD tmp;
846     tmp[0]= valeur[2]*v2[3]-valeur[3]*v2[2];
847     tmp[0]= valeur[3]*v2[0]-valeur[0]*v2[3];
848     tmp[0]= valeur[0]*v2[1]-valeur[1]*v2[0];
849     tmp[0]= valeur[1]*v2[2]-valeur[2]*v2[1];
850     return tmp;
851     }
852 souaissa 69 double2 OT_VECTEUR_4DD::operator* (const OT_VECTEUR_4DD & a)
853     {
854     return a[0]*valeur[0]+a[1]*valeur[1]+a[2]*valeur[2]+a[3]*valeur[3];
855     }
856    
857    
858     OT_VECTEUR_4DD& OT_VECTEUR_4DD::operator+= (const OT_VECTEUR_4DD& rkV)
859     {
860     valeur[0] = valeur[0] +rkV.valeur[0];
861     valeur[1] = valeur[1] +rkV.valeur[1];
862     valeur[2] = valeur[2] +rkV.valeur[2];
863     valeur[3] = valeur[3] +rkV.valeur[3];
864     return *this;
865     }
866    
867     OT_VECTEUR_4DD& OT_VECTEUR_4DD::operator-= (const OT_VECTEUR_4DD& rkV)
868     {
869     valeur[0] = valeur[0] -rkV.valeur[0];
870     valeur[1] = valeur[1] -rkV.valeur[1];
871     valeur[2] = valeur[2] -rkV.valeur[2];
872     valeur[3] = valeur[3] -rkV.valeur[3];
873     return *this;
874     }
875     OT_VECTEUR_4DD& OT_VECTEUR_4DD::operator*= (double2 fScalar)
876     {
877     valeur[0] =valeur[0] * fScalar;
878     valeur[1] =valeur[1] * fScalar;
879     valeur[2] =valeur[2] * fScalar;
880     valeur[3] =valeur[3] * fScalar;
881     return *this;
882     }
883     OT_VECTEUR_4DD& OT_VECTEUR_4DD::operator/= (double2 fScalar)
884     {
885     double2 a = 1/fScalar;
886     valeur[0] =valeur[0] * a;
887     valeur[1] =valeur[1] * a;
888     valeur[2] =valeur[2] * a;
889     valeur[3] =valeur[3] * a;
890     return *this;
891     }
892    
893    
894     std::ostream & operator << ( std::ostream & __os, const OT_VECTEUR_4DD & __vec)
895     {
896 souaissa 70 __os << "[ "<<setw(20)<< __vec[0]<<setw(4)<<", "<<setw(20)<<__vec[1]<<", "<<setw(20)<<__vec[2]<<setw(4)<<", "<<setw(20)<<__vec[3]<<setw(4)<<"] ";
897 souaissa 69 return __os;
898     }
899    
900     int OT_VECTEUR_4DD::compare_valeur (const OT_VECTEUR_4DD& mdd) const
901     {
902 souaissa 75 if (valeur[0]==mdd.valeur[0]&&valeur[1]==mdd.valeur[1] &&valeur[2]==mdd.valeur[2]&&valeur[3]==mdd.valeur[3])
903     return 1;
904     else return 0;
905    
906 souaissa 69 }
907    
908 souaissa 75 bool OT_VECTEUR_4DD::operator==(const OT_VECTEUR_4DD& mdd) const
909 souaissa 69 {
910 francois 76 return compare_valeur(mdd);
911 souaissa 69 }
912    
913    
914    
915    
916    
917    
918    
919    
920    
921    
922    
923    
924    
925    
926    
927    
928    
929    
930    
931    
932 souaissa 58 int OPERATEUR::egal(double a,double b,double eps)
933     {
934     double res=fabs(a-b);
935     if (res<eps) return(1);
936     return(0);
937     }
938    
939    
940     double OPERATEUR::qualite_triangle(double* noeud1,double* noeud2,double* noeud3)
941     {
942     OT_VECTEUR_3D vec12(noeud1,noeud2);
943     OT_VECTEUR_3D vec13(noeud1,noeud3);
944     OT_VECTEUR_3D vec23(noeud2,noeud3);
945     double d12=vec12.get_longueur();
946     double d13=vec13.get_longueur();
947     double d23=vec23.get_longueur();
948     double dmax=max(d12,d13);
949     dmax=max(dmax,d23);
950     double p=0.5*(d12+d13+d23);
951     double crit=(p-d12)*(p-d13)*(p-d23)/p;
952     if (crit<0) crit=0.;
953     crit=sqrt(crit);
954     crit=crit/dmax;
955     crit=crit*3.4641101614;
956     return crit;
957     }
958    
959     void OPERATEUR::doubleto2int(double val,int& val1,int& val2)
960     {
961     int *p=(int*)(&(val));
962     val1=(*p);
963     val2=(*(p+1));
964     }
965    
966     double OPERATEUR::qualite_tetra(double* noeud1,double* noeud2,double* noeud3,double *noeud4)
967     {
968     OT_VECTEUR_3D ab(noeud1,noeud2);
969     OT_VECTEUR_3D ac(noeud1,noeud3);
970     OT_VECTEUR_3D cb(noeud3,noeud2);
971     OT_VECTEUR_3D db(noeud4,noeud2);
972     OT_VECTEUR_3D da(noeud4,noeud1);
973     OT_VECTEUR_3D dc(noeud4,noeud3);
974     double dab=ab.get_longueur();
975     double dac=ac.get_longueur();
976     double dcb=cb.get_longueur();
977     double ddb=db.get_longueur();
978     double dda=da.get_longueur();
979     double ddc=dc.get_longueur();
980    
981     double dmax=max(dab,dac);
982     dmax=max(dmax,dcb);
983     dmax=max(dmax,ddb);
984     dmax=max(dmax,dda);
985     dmax=max(dmax,ddc);
986    
987     OT_VECTEUR_3D pvec=ab&ac;
988     double vol=pvec*da;
989     if (vol>0.) return 0.;
990     vol=fabs(vol);
991     double som=pvec.get_longueur()/2.;
992     pvec=ab&da;
993     som=som+pvec.get_longueur()/2.;
994     pvec=ac&dc;
995     som=som+pvec.get_longueur()/2.;
996     pvec=cb&db;
997     som=som+pvec.get_longueur()/2.;
998     double crit=vol/som/dmax/0.408249;
999     return crit;
1000     }
1001    
1002    
1003    
1004     template <int N> DOUBLEN<N>::DOUBLEN()
1005     {
1006     tab[0]=-1;
1007     }
1008    
1009     template <int N> DOUBLEN<N>::DOUBLEN(DOUBLEN& mdd)
1010     {
1011     for (int i=0;i<N;i++)
1012     tab[i]=mdd.tab[i];
1013     }
1014    
1015     template <int N> DOUBLEN<N>::DOUBLEN(double *v)
1016     {
1017     for (int i=0;i<N;i++)
1018     tab[i]=v[i];
1019     }
1020    
1021     template <int N> DOUBLEN<N>::~DOUBLEN()
1022     {
1023     }
1024    
1025     template <int N> double DOUBLEN<N>::get_valeur(int num)
1026     {
1027     return tab[num];
1028     }
1029    
1030     template <int N> void DOUBLEN<N>::change_valeur(int num,double val)
1031     {
1032     tab[num]=val;
1033     }
1034    
1035     template <int N> DOUBLEN<N> & DOUBLEN<N>::operator=(DOUBLEN<N> & mdd)
1036     {
1037     for (int i=0;i<N;i++)
1038     tab[i]=mdd.tab[i];
1039     return *this;
1040     }
1041    
1042     template <int N> DOUBLEN<N> & DOUBLEN<N>::operator=(const DOUBLEN<N> & mdd)
1043     {
1044     for (int i=0;i<N;i++)
1045     tab[i]=mdd.tab[i];
1046     return *this;
1047     }
1048    
1049    
1050     #ifdef BORLANDCPP
1051     // instanciation pour Windows (Borland c++)
1052     template class __export DOUBLEN<1>;
1053     template class __export DOUBLEN<4>;
1054     template class __export DOUBLEN<10>;
1055     #else
1056     #ifdef GCC
1057     // instanciation pour Unix (g++)
1058     template class DOUBLEN<1>;
1059     template class DOUBLEN<4>;
1060     template class DOUBLEN<10>;
1061     #else
1062     // instanciation pour Windows (Visual Studio)
1063     #endif
1064     #endif
1065    
1066