ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/outil/src/ot_mathematique.cpp
Revision: 19
Committed: Wed Jun 20 17:47:45 2007 UTC (17 years, 10 months ago) by francois
Original Path: magic/lib/outil/outil/src/ot_mathematique.cpp
File size: 20637 byte(s)
Log Message:
ajout d'un index dans les param NURBS

File Contents

# User Rev Content
1 5 //------------------------------------------------------------
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    
27     #include <math.h>
28     #include <algorithm>
29     #include "ot_mathematique.h"
30    
31     using namespace std;
32    
33    
34     OT_VECTEUR_3D::OT_VECTEUR_3D(double x,double y,double z)
35     {
36     valeur[0]=x;
37     valeur[1]=y;
38     valeur[2]=z;
39     }
40    
41    
42     OT_VECTEUR_3D::OT_VECTEUR_3D(double *xyz)
43     {
44     valeur[0]=xyz[0];
45     valeur[1]=xyz[1];
46     valeur[2]=xyz[2];
47     }
48    
49     OT_VECTEUR_3D::OT_VECTEUR_3D(double *xyz1,double *xyz2)
50     {
51     valeur[0]=xyz2[0]-xyz1[0];
52     valeur[1]=xyz2[1]-xyz1[1];
53     valeur[2]=xyz2[2]-xyz1[2];
54     }
55    
56     OT_VECTEUR_3D::OT_VECTEUR_3D(const OT_VECTEUR_3D& mdd)
57     {
58     valeur[0]=mdd.valeur[0];
59     valeur[1]=mdd.valeur[1];
60     valeur[2]=mdd.valeur[2];
61     }
62    
63     OT_VECTEUR_3D::OT_VECTEUR_3D(OT_VECTEUR_3D& mdd)
64     {
65     valeur[0]=mdd.valeur[0];
66     valeur[1]=mdd.valeur[1];
67     valeur[2]=mdd.valeur[2];
68     }
69    
70    
71     double OT_VECTEUR_3D::get_x(void) const
72     {
73     return valeur[0];
74     }
75    
76    
77     double OT_VECTEUR_3D::get_y(void) const
78     {
79     return valeur[1];
80     }
81    
82    
83     double OT_VECTEUR_3D::get_z(void) const
84     {
85     return valeur[2];
86     }
87    
88     void OT_VECTEUR_3D::change_x(double x)
89     {
90     valeur[0]=x;
91     }
92    
93     void OT_VECTEUR_3D::change_y(double y)
94     {
95     valeur[1]=y;
96     }
97    
98     void OT_VECTEUR_3D::change_z(double z)
99     {
100     valeur[2]=z;
101     }
102    
103     double* OT_VECTEUR_3D::get_xyz(void)
104     {
105     return valeur;
106     }
107    
108     OT_VECTEUR_3D::operator const double* () const
109     {
110     return valeur;
111     }
112    
113     OT_VECTEUR_3D::operator double* ()
114     {
115     return valeur;
116     }
117    
118     double OT_VECTEUR_3D::operator[] (int i) const
119     {
120     return valeur[i];
121     }
122    
123     double & OT_VECTEUR_3D::operator[] (int i)
124     {
125     return valeur[i];
126     }
127    
128     double OT_VECTEUR_3D::operator() (int i) const
129     {
130     return valeur[i];
131     }
132     double & OT_VECTEUR_3D::operator() (int i)
133     {
134     return valeur[i];
135     }
136    
137     double OT_VECTEUR_3D::get_longueur(void) const
138     {
139     return sqrt(valeur[0]*valeur[0]+valeur[1]*valeur[1]+valeur[2]*valeur[2]);
140     }
141    
142     double OT_VECTEUR_3D::get_longueur2(void) const
143     {
144     return (valeur[0]*valeur[0]+valeur[1]*valeur[1]+valeur[2]*valeur[2]);
145     }
146    
147    
148     void OT_VECTEUR_3D::norme(void)
149     {
150     double longueur=get_longueur();
151     valeur[0]=valeur[0]/longueur;
152     valeur[2]=valeur[2]/longueur;
153     valeur[1]=valeur[1]/longueur;
154     }
155    
156     int OT_VECTEUR_3D::compare_valeur (const OT_VECTEUR_3D& mdd) const
157     {
158     return memcmp(valeur, mdd.valeur, 3*sizeof(double));
159     }
160    
161     bool OT_VECTEUR_3D::operator== (const OT_VECTEUR_3D& mdd) const
162     {
163     return compare_valeur(mdd) == 0;
164     }
165    
166     bool OT_VECTEUR_3D::operator!= (const OT_VECTEUR_3D& mdd) const
167     {
168     return compare_valeur(mdd) != 0;
169     }
170    
171     bool OT_VECTEUR_3D::operator< (const OT_VECTEUR_3D& mdd) const
172     {
173     return compare_valeur(mdd) < 0;
174     }
175    
176     bool OT_VECTEUR_3D::operator<= (const OT_VECTEUR_3D& mdd) const
177     {
178     return compare_valeur(mdd) <= 0;
179     }
180    
181     bool OT_VECTEUR_3D::operator> (const OT_VECTEUR_3D& mdd) const
182     {
183     return compare_valeur(mdd) > 0;
184     }
185    
186     bool OT_VECTEUR_3D::operator>= (const OT_VECTEUR_3D& mdd) const
187     {
188     return compare_valeur(mdd) >= 0;
189     }
190    
191     std::ostream & operator << ( std::ostream & __os, const OT_VECTEUR_3D & __vec)
192     {
193     __os << "{ "<< __vec[0]<<", "<<__vec[1]<<", "<<__vec[2]<<"}, ";
194     return __os;
195     }
196    
197     OT_VECTEUR_3D operator/ (const OT_VECTEUR_3D& mdd1, double diviseur)
198     {
199     OT_VECTEUR_3D kQuot;
200    
201     if ( diviseur != (double)0.0 )
202     {
203     double facteur = ((double)1.0)/diviseur;
204     kQuot.valeur[0] = facteur*mdd1.valeur[0];
205     kQuot.valeur[1] = facteur*mdd1.valeur[1];
206     kQuot.valeur[2] = facteur*mdd1.valeur[2];
207     }
208     else
209     {
210     kQuot.valeur[0] = 1E300;
211     kQuot.valeur[1] = 1E300;
212     kQuot.valeur[2] = 1E300;
213     }
214    
215     return kQuot;
216     }
217    
218     OT_VECTEUR_3D operator- (const OT_VECTEUR_3D& mdd1)
219     {
220     return OT_VECTEUR_3D(
221     -mdd1.valeur[0],
222     -mdd1.valeur[1],
223     -mdd1.valeur[2]);
224     }
225    
226     double operator*(const OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
227     {
228     double tmp;
229     tmp=mdd1.valeur[0]*mdd2.valeur[0]+mdd1.valeur[1]*mdd2.valeur[1]+mdd1.valeur[2]*mdd2.valeur[2];
230     return tmp;
231     }
232    
233     OT_VECTEUR_3D operator*(double mdd1,const OT_VECTEUR_3D& mdd2)
234     {
235     OT_VECTEUR_3D tmp;
236     tmp.valeur[0]=mdd1*mdd2.valeur[0];
237     tmp.valeur[1]=mdd1*mdd2.valeur[1];
238     tmp.valeur[2]=mdd1*mdd2.valeur[2];
239     return tmp;
240     }
241    
242     OT_VECTEUR_3D& OT_VECTEUR_3D::operator+= (const OT_VECTEUR_3D& mdd)
243     {
244     valeur[0] += mdd.valeur[0];
245     valeur[1] += mdd.valeur[1];
246     valeur[2] += mdd.valeur[2];
247     return *this;
248     }
249     OT_VECTEUR_3D& OT_VECTEUR_3D::operator-= (const OT_VECTEUR_3D& mdd)
250     {
251     valeur[0] -= mdd.valeur[0];
252     valeur[1] -= mdd.valeur[1];
253     valeur[2] -= mdd.valeur[2];
254     return *this;
255     }
256     OT_VECTEUR_3D& OT_VECTEUR_3D::operator*= (double facteur)
257     {
258     valeur[0] *= facteur;
259     valeur[1] *= facteur;
260     valeur[2] *= facteur;
261     return *this;
262     }
263     OT_VECTEUR_3D& OT_VECTEUR_3D::operator/= (double diviseur)
264     {
265     if ( diviseur != (double)0.0 )
266     {
267     double facteur = ((double)1.0)/diviseur;
268     valeur[0] *= facteur;
269     valeur[1] *= facteur;
270     valeur[2] *= facteur;
271     }
272     else
273     {
274     valeur[0] = 1E300;
275     valeur[1] = 1E300;
276     valeur[2] = 1E300;
277     }
278    
279     return *this;
280     }
281    
282     OT_VECTEUR_3D operator&(const OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
283     {
284     OT_VECTEUR_3D tmp;
285     tmp.valeur[0]=mdd1.valeur[1]*mdd2.valeur[2]-mdd1.valeur[2]*mdd2.valeur[1];
286     tmp.valeur[1]=mdd1.valeur[2]*mdd2.valeur[0]-mdd1.valeur[0]*mdd2.valeur[2];
287     tmp.valeur[2]=mdd1.valeur[0]*mdd2.valeur[1]-mdd1.valeur[1]*mdd2.valeur[0];
288     return tmp;
289     }
290    
291     double OT_VECTEUR_3D::diff(void)
292     {
293     return 0.33333333333333*(fabs(valeur[0])+fabs(valeur[1])+fabs(valeur[2]));
294     }
295    
296    
297     OT_MATRICE_3D::OT_MATRICE_3D(OT_VECTEUR_3D& colonne1,OT_VECTEUR_3D& colonne2,OT_VECTEUR_3D& colonne3)
298     {
299     vec[0]=colonne1;
300     vec[1]=colonne2;
301     vec[2]=colonne3;
302     }
303    
304     OT_MATRICE_3D::OT_MATRICE_3D(const OT_MATRICE_3D& mdd)
305     {
306     vec[0]=mdd.vec[0];
307     vec[1]=mdd.vec[1];
308     vec[2]=mdd.vec[2];
309     }
310    
311 francois 19
312 5 std::ostream & operator << ( std::ostream & __os, const OT_MATRICE_3D & __mat)
313     {
314    
315     __os << "{ ";
316     for (unsigned i=0; i<3; i++)
317     {
318     __os << "[ "<< __mat(i,0)<<" ]\t [ "<<__mat(i,1) <<" ]\t [ "<<__mat(i,2);
319     if (i+1 < 3)
320     __os << std::endl << " ";
321     }
322     __os << " }" << std::endl;
323     return __os;
324     }
325 francois 19
326     OT_VECTEUR_3D & OT_MATRICE_3D::operator [](int i)
327     {
328     return vec[i];
329     }
330 5
331 francois 19 double OT_MATRICE_3D::valeur(int iLigne, int iCol) const
332 5 {
333     return vec[iCol][iLigne];
334     }
335     double & OT_MATRICE_3D::valeur(int iLigne, int iCol)
336     {
337     return vec[iCol][iLigne];
338     }
339     double OT_MATRICE_3D::operator() (int iLigne, int iCol) const
340     {
341     return vec[iCol][iLigne];
342     }
343     double & OT_MATRICE_3D::operator() (int iLigne, int iCol)
344     {
345     return vec[iCol][iLigne];
346     }
347     double OT_MATRICE_3D::get_determinant()
348     {
349     double coFacteur0 = valeur(1,1)*valeur(2,2) - valeur(1,2)*valeur(2,1);
350     double coFacteur1 = valeur(1,2)*valeur(2,0) - valeur(1,0)*valeur(2,2);
351     double coFacteur2 = valeur(1,0)*valeur(2,1) - valeur(1,1)*valeur(2,0);
352     double determinant = valeur(0,0)*coFacteur0 + valeur(0,1)*coFacteur1 + valeur(0,2)*coFacteur2;
353     return determinant;
354     }
355     void OT_MATRICE_3D::identite(void)
356     {
357     for (unsigned i = 0; i<3; i++)
358     for (unsigned j = 0; j < 3; j++)
359     vec[i][j] = (i==j) ? 1.0 : 0.0 ;
360     }
361    
362     void OT_MATRICE_3D::transpose(OT_MATRICE_3D& res) const
363     {
364     for (unsigned i = 0; i<3; i++)
365     for (unsigned j = 0; j < 3; j++)
366     res.vec[i][j] = vec[j][i];
367     }
368    
369     OT_MATRICE_3D OT_MATRICE_3D::inverse () const
370     {
371     // Invert a 3x3 using cofactors. This is faster than using a generic
372     // Gaussian elimination because of the loop overhead of such a method.
373    
374     OT_MATRICE_3D matInverse;
375    
376     matInverse.valeur(0,0) = valeur(1,1)*valeur(2,2) - valeur(1,2)*valeur(2,1);
377     matInverse.valeur(0,1) = valeur(0,2)*valeur(2,1) - valeur(0,1)*valeur(2,2);
378     matInverse.valeur(0,2) = valeur(0,1)*valeur(1,2) - valeur(0,2)*valeur(1,1);
379     matInverse.valeur(1,0) = valeur(1,2)*valeur(2,0) - valeur(1,0)*valeur(2,2);
380     matInverse.valeur(1,1) = valeur(0,0)*valeur(2,2) - valeur(0,2)*valeur(2,0);
381     matInverse.valeur(1,2) = valeur(0,2)*valeur(1,0) - valeur(0,0)*valeur(1,2);
382     matInverse.valeur(2,0) = valeur(1,0)*valeur(2,1) - valeur(1,1)*valeur(2,0);
383     matInverse.valeur(2,1) = valeur(0,1)*valeur(2,0) - valeur(0,0)*valeur(2,1);
384     matInverse.valeur(2,2) = valeur(0,0)*valeur(1,1) - valeur(0,1)*valeur(1,0);
385    
386     double det = valeur(0,0)*matInverse(0,0) + valeur(0,1)*matInverse(1,0)+
387     valeur(0,2)*matInverse(2,0);
388    
389     if ( fabs(det) <= 1E-12 )
390     {
391     matInverse.identite();
392     return matInverse;
393     }
394    
395     double invDet = 1/det;
396    
397     for (int i=0; i<3; i++)
398     for (int j=0; j<3; j++)
399     matInverse.vec[i][j] *= invDet;
400    
401     return matInverse;
402     }
403    
404     void OT_MATRICE_3D::change_vecteur1(OT_VECTEUR_3D& v)
405     {
406     vec[0]=v;
407     }
408     void OT_MATRICE_3D::change_vecteur2(OT_VECTEUR_3D& v)
409     {
410     vec[1]=v;
411     }
412     void OT_MATRICE_3D::change_vecteur3(OT_VECTEUR_3D& v)
413     {
414     vec[2]=v;
415     }
416     OT_VECTEUR_3D& OT_MATRICE_3D::get_vecteur1(void)
417     {
418     return vec[0];
419     }
420     OT_VECTEUR_3D& OT_MATRICE_3D::get_vecteur2(void)
421     {
422     return vec[1];
423     }
424     OT_VECTEUR_3D& OT_MATRICE_3D::get_vecteur3(void)
425     {
426     return vec[2];
427     }
428    
429    
430    
431 francois 19 OT_MATRICE_3D OT_MATRICE_3D::VecteursPropres(OT_VECTEUR_3D& Val_Propre,double Eps,int ItMax)
432     {
433     OT_MATRICE_3D m=*this;
434     OT_VECTEUR_3D ZERO(0.,0.,0.) ;
435     OT_MATRICE_3D vct_prpres(ZERO,ZERO,ZERO);
436     OT_VECTEUR_3D V ;
437    
438     int Dim=3;
439     int I,J,K;
440     OT_VECTEUR_3D U, Up;
441     double N,T;
442     //-----
443     for(int vi=0;vi<3;vi++)
444     {
445     //===================
446     Up=vec[0]; // il faut l'initialiser d'une autre maniere
447     K=0;
448     do{
449     U=Up;
450     for(int I=0 ;I<Dim;I++){
451     T=0;
452     for(int J=0; J<Dim;J++) T=T+valeur(I,J)*U[J];
453     Up[I]=T ;
454     }
455     N=Up.get_longueur();
456     if( Up[0]<0) N=-N;
457     T=0;
458     for(int I=0;I<Dim;I++){
459     Up[I]=Up[I]/N;
460     T=T+fabs(Up[I]-U[I]);
461     }
462     K=K+1;
463     }
464     while((K<ItMax)||(T>Eps)) ;
465     V=Up ;
466     vct_prpres[vi]=V;
467     //================
468     double N,D;
469     N=0.0; D=0.0;
470     for(int I=0;I<Dim;I++){
471    
472     D=D+V[I]*V[I];
473     for(int J=0 ;J<Dim;J++) N=N+V[I]*valeur(I,J)*V[J];
474     }
475     double Val=N/D;
476     Val_Propre[vi]=Val;
477     //===============
478     for(int I=0; I<Dim ;I++){
479    
480     for(int J=0;J<Dim;J++) m(I,J)=m(I,J)-Val*V[I]*V[J];
481     }
482    
483     //===============
484    
485     }
486     //-----------------
487    
488    
489     return vct_prpres;
490     }
491    
492     /*
493     double OT_MATRICE_3D::ValeurPropre(double Eps,int ItMax)
494     {
495     double Val_Propre;
496     OT_VECTEUR_3D V=VecteurPropre( Eps, ItMax);
497     int I,J,Dim;
498     double N,D;
499     N=0.0; D=0.0; Dim=3;
500     for(int I=0;I<Dim;I++){
501    
502     D=D+V[I]*V[I];
503     for(int J=0 ;J<Dim;J++) N=N+V[I]*valeur(I,J)*V[J];
504     }
505     Val_Propre=N/D;
506     return Val_Propre;
507     }
508    
509     void OT_MATRICE_3D:: Deflation(OT_VECTEUR_3D& V,double val_prpre)
510     {
511     int Dim=3;
512     for(int I=0; I<Dim ;I++){
513    
514     for(int J=0;J<Dim;J++) valeur(I,J)=valeur(I,J)-val_prpre*V[I]*V[J];
515     }
516    
517     } */
518    
519    
520    
521 5 OT_VECTEUR_3D operator*(const OT_MATRICE_3D& mdd1,const OT_VECTEUR_3D& mdd2)
522     {
523     OT_VECTEUR_3D tmp;
524     tmp(0)=mdd1(0,0)*mdd2(0)+mdd1(0,1)*mdd2(1)+mdd1(0,2)*mdd2(2);
525     tmp(1)=mdd1(1,0)*mdd2(0)+mdd1(1,1)*mdd2(1)+mdd1(1,2)*mdd2(2);
526     tmp(2)=mdd1(2,0)*mdd2(0)+mdd1(2,1)*mdd2(1)+mdd1(2,2)*mdd2(2);
527     return tmp;}
528    
529     OT_VECTEUR_3D operator+(const OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
530     {
531     OT_VECTEUR_3D tmp;
532     tmp(0)=mdd1(0)+mdd2(0);
533     tmp(1)=mdd1(1)+mdd2(1);
534     tmp(2)=mdd1(2)+mdd2(2);
535     return tmp;
536     }
537    
538     OT_VECTEUR_3D operator-(const class OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
539     {
540     OT_VECTEUR_3D tmp;
541     tmp.valeur[0]=mdd1.valeur[0]-mdd2.valeur[0];
542     tmp.valeur[1]=mdd1.valeur[1]-mdd2.valeur[1];
543     tmp.valeur[2]=mdd1.valeur[2]-mdd2.valeur[2];
544     return tmp;
545     }
546    
547     OT_MATRICE_3D operator*(const OT_MATRICE_3D& mdd1,const OT_MATRICE_3D& mdd2)
548     {
549     OT_MATRICE_3D tmp;
550     for (int i=0; i<3; i++)
551     for (int j=0; j<3; j++)
552     {
553     tmp(i,j) = 0;
554     for (int k=0; k<3; k++)
555     tmp(i,j) += mdd1(i,k)*mdd2(k,j);
556     }
557     return tmp;
558     }
559    
560     #ifdef ot_vecteur_4d
561     OT_VECTEUR_4D::OT_VECTEUR_4D (double v[4])
562     {
563     for (int i=0; i<4; i++)
564     valeur[i]=v[i];
565     }
566    
567     OT_VECTEUR_4D::OT_VECTEUR_4D (double x, double y, double z, double w)
568     {
569     valeur[0] = x;
570     valeur[1] = y;
571     valeur[2] = z;
572     valeur[3] = w;
573     }
574    
575     OT_VECTEUR_4D::OT_VECTEUR_4D (const OT_VECTEUR_4D & v)
576     {
577     for (int i=0; i<4; i++)
578     valeur[i]=v[i];
579     }
580    
581    
582     double & OT_VECTEUR_4D::operator()(int i)
583     {
584     return valeur[i];
585     }
586    
587     double OT_VECTEUR_4D::operator()(int i) const
588     {
589     return valeur[i];
590     }
591    
592     double & OT_VECTEUR_4D::operator[](int i)
593     {
594     return valeur[i];
595     }
596    
597     double OT_VECTEUR_4D::operator[](int i) const
598     {
599     return valeur[i];
600     }
601    
602     OT_VECTEUR_4D::operator const double* () const
603     {
604     return valeur;
605     }
606    
607     OT_VECTEUR_4D::operator double* ()
608     {
609     return valeur;
610     }
611    
612     double OT_VECTEUR_4D::get_x(void) const
613     { return valeur[0]; }
614     double & OT_VECTEUR_4D::x(void)
615     { return valeur[0]; }
616     double OT_VECTEUR_4D::get_y(void) const
617     { return valeur[1]; }
618     double & OT_VECTEUR_4D::y(void)
619     { return valeur[1]; }
620     double OT_VECTEUR_4D::get_z(void) const
621     { return valeur[2]; }
622     double & OT_VECTEUR_4D::z(void)
623     { return valeur[2]; }
624     double OT_VECTEUR_4D::get_w(void) const
625     { return valeur[3]; }
626     double & OT_VECTEUR_4D::w(void)
627     { return valeur[3]; }
628     void OT_VECTEUR_4D::change_x(double x)
629     { valeur[0] = x; }
630     void OT_VECTEUR_4D::change_y(double y)
631     { valeur[1] = y; }
632     void OT_VECTEUR_4D::change_z(double z)
633     { valeur[2] = z; }
634     void OT_VECTEUR_4D::change_w(double w)
635     { valeur[3] = w; }
636     double* OT_VECTEUR_4D::get_xyzw(void)
637     { return valeur; }
638    
639     OT_VECTEUR_4D OT_VECTEUR_4D::operator+ (const OT_VECTEUR_4D & a)
640     {
641     return OT_VECTEUR_4D(
642     a(0)+valeur[0],a(1)+valeur[1],a(2)+valeur[2],a(3)+valeur[3]
643     );
644     }
645    
646     OT_VECTEUR_4D OT_VECTEUR_4D::operator- (const OT_VECTEUR_4D & a)
647     {
648     return OT_VECTEUR_4D(
649     valeur[0]-a(0),valeur[1]-a(1),valeur[2]-a(2),valeur[3]-a(3)
650     );
651     }
652    
653     OT_VECTEUR_4D operator* (const OT_VECTEUR_4D & rkV, const double a)
654     {
655     return OT_VECTEUR_4D(
656     a*rkV[0],a*rkV[1],a*rkV[2],a*rkV[3]
657     );
658     }
659     OT_VECTEUR_4D operator* (const double a, const OT_VECTEUR_4D & rkV)
660     {
661     return OT_VECTEUR_4D(
662     a*rkV[0],a*rkV[1],a*rkV[2],a*rkV[3]
663     );
664     }
665    
666     OT_VECTEUR_4D operator/ (const OT_VECTEUR_4D & rkV, const double a)
667     {
668     double inv = 1/a;
669     return OT_VECTEUR_4D(
670     inv*rkV[0],inv*rkV[1],inv*rkV[2],inv*rkV[3]
671     );
672     }
673    
674    
675     OT_VECTEUR_4D operator/ (const double a, const OT_VECTEUR_4D & rkV)
676     {
677     double inv = 1/a;
678     return OT_VECTEUR_4D(
679     inv*rkV[0],inv*rkV[1],inv*rkV[2],inv*rkV[3]
680     );
681     }
682    
683     double OT_VECTEUR_4D::operator* (const OT_VECTEUR_4D & a)
684     {
685     return a[0]*valeur[0]+a[1]*valeur[1]+a[2]*valeur[2]+a[3]*valeur[3];
686     }
687    
688    
689     OT_VECTEUR_4D& OT_VECTEUR_4D::operator+= (const OT_VECTEUR_4D& rkV)
690     {
691     valeur[0] += rkV.valeur[0];
692     valeur[1] += rkV.valeur[1];
693     valeur[2] += rkV.valeur[2];
694     valeur[3] += rkV.valeur[3];
695     return *this;
696     }
697    
698     OT_VECTEUR_4D& OT_VECTEUR_4D::operator-= (const OT_VECTEUR_4D& rkV)
699     {
700     valeur[0] -= rkV.valeur[0];
701     valeur[1] -= rkV.valeur[1];
702     valeur[2] -= rkV.valeur[2];
703     valeur[3] -= rkV.valeur[3];
704     return *this;
705     }
706     OT_VECTEUR_4D& OT_VECTEUR_4D::operator*= (double fScalar)
707     {
708     valeur[0] *= fScalar;
709     valeur[1] *= fScalar;
710     valeur[2] *= fScalar;
711     valeur[3] *= fScalar;
712     return *this;
713     }
714     OT_VECTEUR_4D& OT_VECTEUR_4D::operator/= (double fScalar)
715     {
716     double a = 1/fScalar;
717     valeur[0] *= a;
718     valeur[1] *= a;
719     valeur[2] *= a;
720     valeur[3] *= a;
721     return *this;
722     }
723     OT_VECTEUR_4D &OT_VECTEUR_4D::operator= (double a)
724     {
725     valeur[0]=valeur[1]=valeur[2]=valeur[3]=a;
726     return *this;
727     }
728    
729     std::ostream & operator << ( std::ostream & __os, const OT_VECTEUR_4D & __vec)
730     {
731     __os << "{ "<< __vec[0]<<", "<<__vec[1]<<", "<<__vec[2]<<", "<<__vec[3]<<"}, ";
732     return __os;
733     }
734     #endif // ot_vecteur_4d
735    
736     int OPERATEUR::egal(double a,double b,double eps)
737     {
738     double res=fabs(a-b);
739     if (res<eps) return(1);
740     return(0);
741     }
742    
743    
744     double OPERATEUR::qualite_triangle(double* noeud1,double* noeud2,double* noeud3)
745     {
746     OT_VECTEUR_3D vec12(noeud1,noeud2);
747     OT_VECTEUR_3D vec13(noeud1,noeud3);
748     OT_VECTEUR_3D vec23(noeud2,noeud3);
749     double d12=vec12.get_longueur();
750     double d13=vec13.get_longueur();
751     double d23=vec23.get_longueur();
752     double dmax=max(d12,d13);
753     dmax=max(dmax,d23);
754     double p=0.5*(d12+d13+d23);
755     double crit=(p-d12)*(p-d13)*(p-d23)/p;
756     if (crit<0) crit=0.;
757     crit=sqrt(crit);
758     crit=crit/dmax;
759     crit=crit*3.4641101614;
760     return crit;
761     }
762    
763     void OPERATEUR::doubleto2int(double val,int& val1,int& val2)
764     {
765     int *p=(int*)(&(val));
766     val1=(*p);
767     val2=(*(p+1));
768     }
769    
770     double OPERATEUR::qualite_tetra(double* noeud1,double* noeud2,double* noeud3,double *noeud4)
771     {
772     OT_VECTEUR_3D ab(noeud1,noeud2);
773     OT_VECTEUR_3D ac(noeud1,noeud3);
774     OT_VECTEUR_3D cb(noeud3,noeud2);
775     OT_VECTEUR_3D db(noeud4,noeud2);
776     OT_VECTEUR_3D da(noeud4,noeud1);
777     OT_VECTEUR_3D dc(noeud4,noeud3);
778     double dab=ab.get_longueur();
779     double dac=ac.get_longueur();
780     double dcb=cb.get_longueur();
781     double ddb=db.get_longueur();
782     double dda=da.get_longueur();
783     double ddc=dc.get_longueur();
784    
785     double dmax=max(dab,dac);
786     dmax=max(dmax,dcb);
787     dmax=max(dmax,ddb);
788     dmax=max(dmax,dda);
789     dmax=max(dmax,ddc);
790    
791     OT_VECTEUR_3D pvec=ab&ac;
792     double vol=pvec*da;
793     if (vol>0.) return 0.;
794     vol=fabs(vol);
795     double som=pvec.get_longueur()/2.;
796     pvec=ab&da;
797     som=som+pvec.get_longueur()/2.;
798     pvec=ac&dc;
799     som=som+pvec.get_longueur()/2.;
800     pvec=cb&db;
801     som=som+pvec.get_longueur()/2.;
802     double crit=vol/som/dmax/0.408249;
803     return crit;
804     }
805    
806    
807    
808     template <int N> DOUBLEN<N>::DOUBLEN()
809     {
810     tab[0]=-1;
811     }
812    
813     template <int N> DOUBLEN<N>::DOUBLEN(DOUBLEN& mdd)
814     {
815     for (int i=0;i<N;i++)
816     tab[i]=mdd.tab[i];
817     }
818    
819     template <int N> DOUBLEN<N>::DOUBLEN(double *v)
820     {
821     for (int i=0;i<N;i++)
822     tab[i]=v[i];
823     }
824    
825     template <int N> DOUBLEN<N>::~DOUBLEN()
826     {
827     }
828    
829     template <int N> double DOUBLEN<N>::get_valeur(int num)
830     {
831     return tab[num];
832     }
833    
834     template <int N> void DOUBLEN<N>::change_valeur(int num,double val)
835     {
836     tab[num]=val;
837     }
838    
839     template <int N> DOUBLEN<N> & DOUBLEN<N>::operator=(DOUBLEN<N> & mdd)
840     {
841     for (int i=0;i<N;i++)
842     tab[i]=mdd.tab[i];
843     return *this;
844     }
845    
846     template <int N> DOUBLEN<N> & DOUBLEN<N>::operator=(const DOUBLEN<N> & mdd)
847     {
848     for (int i=0;i<N;i++)
849     tab[i]=mdd.tab[i];
850     return *this;
851     }
852    
853    
854     #ifdef BORLANDCPP
855     // instanciation pour Windows (Borland c++)
856     template class __export DOUBLEN<1>;
857     template class __export DOUBLEN<4>;
858     template class __export DOUBLEN<10>;
859     #else
860     #ifdef GCC
861     // instanciation pour Unix (g++)
862     template class DOUBLEN<1>;
863     template class DOUBLEN<4>;
864     template class DOUBLEN<10>;
865     #else
866     // instanciation pour Windows (Visual Studio)
867     #endif
868     #endif
869    
870