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

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     // 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     #include <iomanip>
27     #include <math.h>
28     #include <string.h>
29     #include <assert.h>
30     #include <algorithm>
31     #include "ot_mathematique.h"
32    
33     using namespace std;
34    
35 couturad 906 OT_QUATERNION::OT_QUATERNION(void)
36     {
37     valeur[0]=0;
38     valeur[1]=0;
39     valeur[2]=0;
40     valeur[3]=0;
41     }
42 francois 283
43 couturad 906 OT_QUATERNION::OT_QUATERNION(double x, double y, double z, double w)
44     {
45     valeur[0]=x;
46     valeur[1]=y;
47     valeur[2]=z;
48     valeur[3]=w;
49     }
50    
51     OT_QUATERNION::OT_QUATERNION(OT_QUATERNION& mdd)
52     {
53     valeur[0]=mdd.valeur[0];
54     valeur[1]=mdd.valeur[1];
55     valeur[2]=mdd.valeur[2];
56     valeur[3]=mdd.valeur[3];
57     }
58    
59     double OT_QUATERNION::get_x(void)
60     {
61     return valeur[0];
62     }
63    
64     void OT_QUATERNION::change_x(double x)
65     {
66     valeur[0]=x;
67     }
68    
69     double OT_QUATERNION::get_y(void)
70     {
71     return valeur[1];
72     }
73    
74     void OT_QUATERNION::change_y(double y)
75     {
76     valeur[1]=y;
77     }
78    
79     double OT_QUATERNION::get_z(void)
80     {
81     return valeur[2];
82     }
83    
84     void OT_QUATERNION::change_z(double z)
85     {
86     valeur[2]=z;
87     }
88    
89     double OT_QUATERNION::get_w(void)
90     {
91     return valeur[3];
92     }
93    
94     void OT_QUATERNION::change_w(double w)
95     {
96     valeur[3]=w;
97     }
98    
99     OT_QUATERNION& OT_QUATERNION::operator=(OT_QUATERNION& mdd)
100     {
101     valeur[0]=mdd.valeur[0];
102     valeur[1]=mdd.valeur[1];
103     valeur[2]=mdd.valeur[2];
104     valeur[3]=mdd.valeur[3];
105     return *this;
106     }
107    
108    
109 francois 283 OT_VECTEUR_3D::OT_VECTEUR_3D(double x,double y,double z)
110 foucault 27 {
111 francois 283 valeur[0]=x;
112     valeur[1]=y;
113     valeur[2]=z;
114 foucault 27 }
115 francois 249
116    
117 francois 283 OT_VECTEUR_3D::OT_VECTEUR_3D(double *xyz)
118     {
119     valeur[0]=xyz[0];
120     valeur[1]=xyz[1];
121     valeur[2]=xyz[2];
122     }
123 francois 249
124 francois 283 OT_VECTEUR_3D::OT_VECTEUR_3D(double *xyz1,double *xyz2)
125     {
126     valeur[0]=xyz2[0]-xyz1[0];
127     valeur[1]=xyz2[1]-xyz1[1];
128     valeur[2]=xyz2[2]-xyz1[2];
129     }
130 souaissa 69
131 francois 283 OT_VECTEUR_3D::OT_VECTEUR_3D(const OT_VECTEUR_3D& mdd)
132     {
133     valeur[0]=mdd.valeur[0];
134     valeur[1]=mdd.valeur[1];
135     valeur[2]=mdd.valeur[2];
136     }
137    
138     OT_VECTEUR_3D::OT_VECTEUR_3D(void)
139     {
140     valeur[0]=0.;
141     valeur[1]=0.;
142     valeur[2]=0.;
143     }
144    
145     double OT_VECTEUR_3D::get_x(void) const
146     {
147     return valeur[0];
148     }
149    
150    
151     double OT_VECTEUR_3D::get_y(void) const
152     {
153     return valeur[1];
154     }
155    
156    
157     double OT_VECTEUR_3D::get_z(void) const
158     {
159     return valeur[2];
160     }
161    
162     void OT_VECTEUR_3D::change_x(double x)
163     {
164     valeur[0]=x;
165     }
166    
167     void OT_VECTEUR_3D::change_y(double y)
168     {
169     valeur[1]=y;
170     }
171    
172     void OT_VECTEUR_3D::change_z(double z)
173     {
174     valeur[2]=z;
175     }
176    
177     double* OT_VECTEUR_3D::get_xyz(void)
178     {
179     return valeur;
180     }
181    
182     OT_VECTEUR_3D::operator const double* () const
183     {
184     return valeur;
185     }
186    
187     OT_VECTEUR_3D::operator double* ()
188     {
189     return valeur;
190     }
191    
192     double OT_VECTEUR_3D::operator[] (int i) const
193     {
194     return valeur[i];
195     }
196    
197     double & OT_VECTEUR_3D::operator[] (int i)
198     {
199     return valeur[i];
200     }
201    
202     double OT_VECTEUR_3D::operator() (int i) const
203     {
204     return valeur[i];
205     }
206     double & OT_VECTEUR_3D::operator() (int i)
207     {
208     return valeur[i];
209     }
210    
211     double OT_VECTEUR_3D::get_longueur(void) const
212     {
213     return sqrt(valeur[0]*valeur[0]+valeur[1]*valeur[1]+valeur[2]*valeur[2]);
214     }
215    
216     double OT_VECTEUR_3D::get_longueur2(void) const
217     {
218     return (valeur[0]*valeur[0]+valeur[1]*valeur[1]+valeur[2]*valeur[2]);
219     }
220    
221    
222     void OT_VECTEUR_3D::norme(void)
223     {
224     double longueur=get_longueur();
225     if (longueur != 0)
226     {
227     valeur[0]=valeur[0]/longueur;
228     valeur[2]=valeur[2]/longueur;
229     valeur[1]=valeur[1]/longueur;
230     }
231     }
232    
233     OT_VECTEUR_3D OT_VECTEUR_3D::gram_shmidt(const OT_VECTEUR_3D& vint)
234     {
235     OT_VECTEUR_3D v(0.,0.,0.);
236     double alpha=(*this)*vint;
237     v= *this-alpha*vint;
238     return v;
239     }
240    
241     OT_VECTEUR_3D OT_VECTEUR_3D::unite(int i)
242     {
243     OT_VECTEUR_3D v(0.,0.,0.);
244     v[i]=1.;
245     return v;
246     }
247    
248    
249     int OT_VECTEUR_3D::compare_valeur (const OT_VECTEUR_3D& mdd) const
250     {
251     return memcmp(valeur, mdd.valeur, 3*sizeof(double));
252     }
253    
254     bool OT_VECTEUR_3D::operator== (const OT_VECTEUR_3D& mdd) const
255     {
256     return compare_valeur(mdd) == 0;
257     }
258    
259     bool OT_VECTEUR_3D::operator!= (const OT_VECTEUR_3D& mdd) const
260     {
261     return compare_valeur(mdd) != 0;
262     }
263    
264     bool OT_VECTEUR_3D::operator< (const OT_VECTEUR_3D& mdd) const
265     {
266     return compare_valeur(mdd) < 0;
267     }
268    
269     bool OT_VECTEUR_3D::operator<= (const OT_VECTEUR_3D& mdd) const
270     {
271     return compare_valeur(mdd) <= 0;
272     }
273    
274     bool OT_VECTEUR_3D::operator> (const OT_VECTEUR_3D& mdd) const
275     {
276     return compare_valeur(mdd) > 0;
277     }
278    
279     bool OT_VECTEUR_3D::operator>= (const OT_VECTEUR_3D& mdd) const
280     {
281     return compare_valeur(mdd) >= 0;
282     }
283    
284     std::ostream & operator << ( std::ostream & __os, const OT_VECTEUR_3D & __vec)
285     {
286     __os << "{ "<< __vec[0]<<", "<<__vec[1]<<", "<<__vec[2]<<"}, ";
287     return __os;
288     }
289    
290     OT_VECTEUR_3D operator/ (const OT_VECTEUR_3D& mdd1, double diviseur)
291     {
292     OT_VECTEUR_3D kQuot;
293    
294     if ( diviseur != (double)0.0 )
295     {
296     double facteur = ((double)1.0)/diviseur;
297     kQuot.valeur[0] = facteur*mdd1.valeur[0];
298     kQuot.valeur[1] = facteur*mdd1.valeur[1];
299     kQuot.valeur[2] = facteur*mdd1.valeur[2];
300     }
301     else
302     {
303     kQuot.valeur[0] = 1E300;
304     kQuot.valeur[1] = 1E300;
305     kQuot.valeur[2] = 1E300;
306     }
307    
308     return kQuot;
309     }
310    
311     OT_VECTEUR_3D operator- (const OT_VECTEUR_3D& mdd1)
312     {
313     return OT_VECTEUR_3D(
314     -mdd1.valeur[0],
315     -mdd1.valeur[1],
316     -mdd1.valeur[2]);
317     }
318    
319     double operator*(const OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
320     {
321     double tmp;
322     tmp=mdd1.valeur[0]*mdd2.valeur[0]+mdd1.valeur[1]*mdd2.valeur[1]+mdd1.valeur[2]*mdd2.valeur[2];
323     return tmp;
324     }
325    
326     OT_VECTEUR_3D operator*(double mdd1,const OT_VECTEUR_3D& mdd2)
327     {
328     OT_VECTEUR_3D tmp;
329     tmp.valeur[0]=mdd1*mdd2.valeur[0];
330     tmp.valeur[1]=mdd1*mdd2.valeur[1];
331     tmp.valeur[2]=mdd1*mdd2.valeur[2];
332     return tmp;
333     }
334    
335     OT_VECTEUR_3D& OT_VECTEUR_3D::operator+= (const OT_VECTEUR_3D& mdd)
336     {
337     valeur[0] += mdd.valeur[0];
338     valeur[1] += mdd.valeur[1];
339     valeur[2] += mdd.valeur[2];
340     return *this;
341     }
342     OT_VECTEUR_3D& OT_VECTEUR_3D::operator-= (const OT_VECTEUR_3D& mdd)
343     {
344     valeur[0] -= mdd.valeur[0];
345     valeur[1] -= mdd.valeur[1];
346     valeur[2] -= mdd.valeur[2];
347     return *this;
348     }
349     OT_VECTEUR_3D& OT_VECTEUR_3D::operator*= (double facteur)
350     {
351     valeur[0] *= facteur;
352     valeur[1] *= facteur;
353     valeur[2] *= facteur;
354     return *this;
355     }
356     OT_VECTEUR_3D& OT_VECTEUR_3D::operator/= (double diviseur)
357     {
358     if ( diviseur != (double)0.0 )
359     {
360     double facteur = ((double)1.0)/diviseur;
361     valeur[0] *= facteur;
362     valeur[1] *= facteur;
363     valeur[2] *= facteur;
364     }
365     else
366     {
367     valeur[0] = 1E300;
368     valeur[1] = 1E300;
369     valeur[2] = 1E300;
370     }
371    
372     return *this;
373     }
374    
375     OT_VECTEUR_3D operator&(const OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
376     {
377     OT_VECTEUR_3D tmp;
378     tmp.valeur[0]=mdd1.valeur[1]*mdd2.valeur[2]-mdd1.valeur[2]*mdd2.valeur[1];
379     tmp.valeur[1]=mdd1.valeur[2]*mdd2.valeur[0]-mdd1.valeur[0]*mdd2.valeur[2];
380     tmp.valeur[2]=mdd1.valeur[0]*mdd2.valeur[1]-mdd1.valeur[1]*mdd2.valeur[0];
381     return tmp;
382     }
383    
384     double OT_VECTEUR_3D::diff(void)
385     {
386     return 0.33333333333333*(fabs(valeur[0])+fabs(valeur[1])+fabs(valeur[2]));
387     }
388    
389    
390     OT_MATRICE_3D::OT_MATRICE_3D(OT_VECTEUR_3D& colonne1,OT_VECTEUR_3D& colonne2,OT_VECTEUR_3D& colonne3)
391     {
392     vec[0]=colonne1;
393     vec[1]=colonne2;
394     vec[2]=colonne3;
395     }
396    
397     OT_MATRICE_3D::OT_MATRICE_3D(const OT_MATRICE_3D& mdd)
398     {
399     vec[0]=mdd.vec[0];
400     vec[1]=mdd.vec[1];
401     vec[2]=mdd.vec[2];
402     }
403    
404     OT_MATRICE_3D::OT_MATRICE_3D(void)
405     {
406     for (int i=0;i<3;i++)
407     {
408     for (int j=0;j<3;j++)
409     {
410     vec[i][j]=0.;
411     }
412    
413     }
414    
415     }
416     std::ostream & operator << ( std::ostream & __os, const OT_MATRICE_3D & __mat)
417     {
418    
419     __os << "{ ";
420     for (unsigned i=0; i<3; i++)
421     {
422     __os << "[ "<< __mat(i,0)<<" ]\t [ "<<__mat(i,1) <<" ]\t [ "<<__mat(i,2);
423     if (i+1 < 3)
424     __os << std::endl << " ";
425     }
426     __os << " }" << std::endl;
427     return __os;
428     }
429    
430     OT_VECTEUR_3D & OT_MATRICE_3D::operator [](int i)
431     {
432     return vec[i];
433     }
434    
435     OT_MATRICE_3D::OT_MATRICE_3D(double* t)
436     {
437     for (int i=0;i<3;i++)
438     {
439     for (int j=0;j<3;j++)
440     {
441     vec[i][j]=t[i*3+j];
442     }
443    
444     }
445     }
446    
447     double OT_MATRICE_3D::valeur(int iLigne, int iCol) const
448     {
449     return vec[iCol][iLigne];
450     }
451     double & OT_MATRICE_3D::valeur(int iLigne, int iCol)
452     {
453     return vec[iCol][iLigne];
454     }
455     double OT_MATRICE_3D::operator() (int iLigne, int iCol) const
456     {
457     return vec[iCol][iLigne];
458     }
459     double & OT_MATRICE_3D::operator() (int iLigne, int iCol)
460     {
461     return vec[iCol][iLigne];
462     }
463     double OT_MATRICE_3D::get_determinant()
464     {
465     double coFacteur0 = valeur(1,1)*valeur(2,2) - valeur(1,2)*valeur(2,1);
466     double coFacteur1 = valeur(1,2)*valeur(2,0) - valeur(1,0)*valeur(2,2);
467     double coFacteur2 = valeur(1,0)*valeur(2,1) - valeur(1,1)*valeur(2,0);
468     double determinant = valeur(0,0)*coFacteur0 + valeur(0,1)*coFacteur1 + valeur(0,2)*coFacteur2;
469     return determinant;
470     }
471     void OT_MATRICE_3D::identite(void)
472     {
473     for (unsigned i = 0; i<3; i++)
474     for (unsigned j = 0; j < 3; j++)
475     vec[i][j] = (i==j) ? 1.0 : 0.0 ;
476     }
477    
478     void OT_MATRICE_3D::transpose(OT_MATRICE_3D& res) const
479     {
480     for (unsigned i = 0; i<3; i++)
481     for (unsigned j = 0; j < 3; j++)
482     res.vec[i][j] = vec[j][i];
483     }
484    
485     OT_MATRICE_3D OT_MATRICE_3D::transpose() const
486     {
487     OT_MATRICE_3D temp;
488     for (unsigned i = 0; i<3; i++)
489     for (unsigned j = 0; j < 3; j++)
490     temp.valeur(i,j) =valeur(j,i);
491    
492     return temp;
493     }
494    
495     OT_MATRICE_3D OT_MATRICE_3D::inverse () const
496     {
497     // Invert a 3x3 using cofactors. This is faster than using a generic
498     // Gaussian elimination because of the loop overhead of such a method.
499    
500     OT_MATRICE_3D matInverse;
501    
502     matInverse.valeur(0,0) = valeur(1,1)*valeur(2,2) - valeur(1,2)*valeur(2,1);
503     matInverse.valeur(0,1) = valeur(0,2)*valeur(2,1) - valeur(0,1)*valeur(2,2);
504     matInverse.valeur(0,2) = valeur(0,1)*valeur(1,2) - valeur(0,2)*valeur(1,1);
505     matInverse.valeur(1,0) = valeur(1,2)*valeur(2,0) - valeur(1,0)*valeur(2,2);
506     matInverse.valeur(1,1) = valeur(0,0)*valeur(2,2) - valeur(0,2)*valeur(2,0);
507     matInverse.valeur(1,2) = valeur(0,2)*valeur(1,0) - valeur(0,0)*valeur(1,2);
508     matInverse.valeur(2,0) = valeur(1,0)*valeur(2,1) - valeur(1,1)*valeur(2,0);
509     matInverse.valeur(2,1) = valeur(0,1)*valeur(2,0) - valeur(0,0)*valeur(2,1);
510     matInverse.valeur(2,2) = valeur(0,0)*valeur(1,1) - valeur(0,1)*valeur(1,0);
511    
512     double det = valeur(0,0)*matInverse(0,0) + valeur(0,1)*matInverse(1,0)+
513     valeur(0,2)*matInverse(2,0);
514    
515     if ( fabs(det) <= 1E-100 )
516     {
517     matInverse.identite();
518     return matInverse;
519     }
520    
521     double invDet = 1/det;
522    
523     for (int i=0; i<3; i++)
524     for (int j=0; j<3; j++)
525     matInverse.vec[i][j] *= invDet;
526    
527     return matInverse;
528     }
529    
530     void OT_MATRICE_3D::change_vecteur1(OT_VECTEUR_3D v)
531     {
532     vec[0]=v;
533     }
534     void OT_MATRICE_3D::change_vecteur2(OT_VECTEUR_3D v)
535     {
536     vec[1]=v;
537     }
538     void OT_MATRICE_3D::change_vecteur3(OT_VECTEUR_3D v)
539     {
540     vec[2]=v;
541     }
542     OT_VECTEUR_3D& OT_MATRICE_3D::get_vecteur1(void)
543     {
544     return vec[0];
545     }
546     OT_VECTEUR_3D& OT_MATRICE_3D::get_vecteur2(void)
547     {
548     return vec[1];
549     }
550     OT_VECTEUR_3D& OT_MATRICE_3D::get_vecteur3(void)
551     {
552     return vec[2];
553     }
554    
555    
556    
557     OT_MATRICE_3D operator+(const OT_MATRICE_3D& mdd1,const OT_MATRICE_3D& mdd2)
558     {
559     OT_MATRICE_3D tmp;
560     for (int i=0;i<3;i++) {
561     for (int j=0;j<3;j++) tmp(i,j)= mdd1(i,j)+mdd2(i,j);
562     }
563     return tmp;
564     }
565    
566    
567     OT_VECTEUR_3D operator*(const OT_MATRICE_3D& mdd1,const OT_VECTEUR_3D& mdd2)
568     {
569     OT_VECTEUR_3D tmp;
570     tmp(0)=mdd1(0,0)*mdd2(0)+mdd1(0,1)*mdd2(1)+mdd1(0,2)*mdd2(2);
571     tmp(1)=mdd1(1,0)*mdd2(0)+mdd1(1,1)*mdd2(1)+mdd1(1,2)*mdd2(2);
572     tmp(2)=mdd1(2,0)*mdd2(0)+mdd1(2,1)*mdd2(1)+mdd1(2,2)*mdd2(2);
573     return tmp;
574     }
575    
576     OT_VECTEUR_3D operator+(const OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
577     {
578     OT_VECTEUR_3D tmp;
579     tmp(0)=mdd1(0)+mdd2(0);
580     tmp(1)=mdd1(1)+mdd2(1);
581     tmp(2)=mdd1(2)+mdd2(2);
582     return tmp;
583     }
584    
585     OT_VECTEUR_3D operator-(const class OT_VECTEUR_3D& mdd1,const OT_VECTEUR_3D& mdd2)
586     {
587     OT_VECTEUR_3D tmp;
588     tmp.valeur[0]=mdd1.valeur[0]-mdd2.valeur[0];
589     tmp.valeur[1]=mdd1.valeur[1]-mdd2.valeur[1];
590     tmp.valeur[2]=mdd1.valeur[2]-mdd2.valeur[2];
591     return tmp;
592     }
593    
594     OT_MATRICE_3D operator*(const OT_MATRICE_3D& mdd1,const OT_MATRICE_3D& mdd2)
595     {
596     OT_MATRICE_3D tmp;
597     for (int i=0; i<3; i++)
598     for (int j=0; j<3; j++)
599     {
600     tmp(i,j) = 0;
601     for (int k=0; k<3; k++)
602     tmp(i,j) += mdd1(i,k)*mdd2(k,j);
603     }
604     return tmp;
605     }
606    
607     #ifdef ot_vecteur_4d
608     OT_VECTEUR_4D::OT_VECTEUR_4D (double v[4])
609     {
610     for (int i=0; i<4; i++)
611     valeur[i]=v[i];
612     }
613    
614     OT_VECTEUR_4D::OT_VECTEUR_4D (double x, double y, double z, double w)
615     {
616     valeur[0] = x;
617     valeur[1] = y;
618     valeur[2] = z;
619     valeur[3] = w;
620     }
621    
622     OT_VECTEUR_4D::OT_VECTEUR_4D (const OT_VECTEUR_4D & v)
623     {
624     for (int i=0; i<4; i++)
625     valeur[i]=v[i];
626     }
627    
628    
629     double & OT_VECTEUR_4D::operator()(int i)
630     {
631     return valeur[i];
632     }
633    
634     double OT_VECTEUR_4D::operator()(int i) const
635     {
636     return valeur[i];
637     }
638    
639     double & OT_VECTEUR_4D::operator[](int i)
640     {
641     return valeur[i];
642     }
643    
644     double OT_VECTEUR_4D::operator[](int i) const
645     {
646     return valeur[i];
647     }
648    
649     OT_VECTEUR_4D::operator const double* () const
650     {
651     return valeur;
652     }
653    
654     OT_VECTEUR_4D::operator double* ()
655     {
656     return valeur;
657     }
658    
659     double OT_VECTEUR_4D::get_x(void) const
660     {
661     return valeur[0];
662     }
663     double & OT_VECTEUR_4D::x(void)
664     {
665     return valeur[0];
666     }
667     double OT_VECTEUR_4D::get_y(void) const
668     {
669     return valeur[1];
670     }
671     double & OT_VECTEUR_4D::y(void)
672     {
673     return valeur[1];
674     }
675     double OT_VECTEUR_4D::get_z(void) const
676     {
677     return valeur[2];
678     }
679     double & OT_VECTEUR_4D::z(void)
680     {
681     return valeur[2];
682     }
683     double OT_VECTEUR_4D::get_w(void) const
684     {
685     return valeur[3];
686     }
687     double & OT_VECTEUR_4D::w(void)
688     {
689     return valeur[3];
690     }
691     void OT_VECTEUR_4D::change_x(double x)
692     {
693     valeur[0] = x;
694     }
695     void OT_VECTEUR_4D::change_y(double y)
696     {
697     valeur[1] = y;
698     }
699     void OT_VECTEUR_4D::change_z(double z)
700     {
701     valeur[2] = z;
702     }
703     void OT_VECTEUR_4D::change_w(double w)
704     {
705     valeur[3] = w;
706     }
707     double* OT_VECTEUR_4D::get_xyzw(void)
708     {
709     return valeur;
710     }
711    
712     OT_VECTEUR_4D OT_VECTEUR_4D::operator+ (const OT_VECTEUR_4D & a)
713     {
714     return OT_VECTEUR_4D(
715     a(0)+valeur[0],a(1)+valeur[1],a(2)+valeur[2],a(3)+valeur[3]
716     );
717     }
718    
719     OT_VECTEUR_4D OT_VECTEUR_4D::operator- (const OT_VECTEUR_4D & a)
720     {
721     return OT_VECTEUR_4D(
722     valeur[0]-a(0),valeur[1]-a(1),valeur[2]-a(2),valeur[3]-a(3)
723     );
724     }
725    
726     OT_VECTEUR_4D operator* (const OT_VECTEUR_4D & rkV, const double a)
727     {
728     return OT_VECTEUR_4D(
729     a*rkV[0],a*rkV[1],a*rkV[2],a*rkV[3]
730     );
731     }
732     OT_VECTEUR_4D operator* (const double a, const OT_VECTEUR_4D & rkV)
733     {
734     return OT_VECTEUR_4D(
735     a*rkV[0],a*rkV[1],a*rkV[2],a*rkV[3]
736     );
737     }
738    
739     OT_VECTEUR_4D operator/ (const OT_VECTEUR_4D & rkV, const double a)
740     {
741     double inv = 1/a;
742     return OT_VECTEUR_4D(
743     inv*rkV[0],inv*rkV[1],inv*rkV[2],inv*rkV[3]
744     );
745     }
746    
747    
748     OT_VECTEUR_4D operator/ (const double a, const OT_VECTEUR_4D & rkV)
749     {
750     double inv = 1/a;
751     return OT_VECTEUR_4D(
752     inv*rkV[0],inv*rkV[1],inv*rkV[2],inv*rkV[3]
753     );
754     }
755    
756     double OT_VECTEUR_4D::operator* (const OT_VECTEUR_4D & a)
757     {
758     return a[0]*valeur[0]+a[1]*valeur[1]+a[2]*valeur[2]+a[3]*valeur[3];
759     }
760    
761    
762     OT_VECTEUR_4D& OT_VECTEUR_4D::operator+= (const OT_VECTEUR_4D& rkV)
763     {
764     valeur[0] += rkV.valeur[0];
765     valeur[1] += rkV.valeur[1];
766     valeur[2] += rkV.valeur[2];
767     valeur[3] += rkV.valeur[3];
768     return *this;
769     }
770    
771     OT_VECTEUR_4D& OT_VECTEUR_4D::operator-= (const OT_VECTEUR_4D& rkV)
772     {
773     valeur[0] -= rkV.valeur[0];
774     valeur[1] -= rkV.valeur[1];
775     valeur[2] -= rkV.valeur[2];
776     valeur[3] -= rkV.valeur[3];
777     return *this;
778     }
779     OT_VECTEUR_4D& OT_VECTEUR_4D::operator*= (double fScalar)
780     {
781     valeur[0] *= fScalar;
782     valeur[1] *= fScalar;
783     valeur[2] *= fScalar;
784     valeur[3] *= fScalar;
785     return *this;
786     }
787     OT_VECTEUR_4D& OT_VECTEUR_4D::operator/= (double fScalar)
788     {
789     double a = 1/fScalar;
790     valeur[0] *= a;
791     valeur[1] *= a;
792     valeur[2] *= a;
793     valeur[3] *= a;
794     return *this;
795     }
796     OT_VECTEUR_4D &OT_VECTEUR_4D::operator= (double a)
797     {
798     valeur[0]=valeur[1]=valeur[2]=valeur[3]=a;
799     return *this;
800     }
801    
802     std::ostream & operator << ( std::ostream & __os, const OT_VECTEUR_4D & __vec)
803     {
804     __os << "{ "<< __vec[0]<<", "<<__vec[1]<<", "<<__vec[2]<<", "<<__vec[3]<<"}, ";
805     return __os;
806     }
807     #endif // ot_vecteur_4d
808    
809    
810    
811     // ot_vecteur_3d
812    
813    
814    
815     OT_VECTEUR_3DD::OT_VECTEUR_3DD (double2 v[3])
816     {
817     for (int i=0; i<3; i++)
818     valeur[i]=v[i];
819     }
820    
821     OT_VECTEUR_3DD::OT_VECTEUR_3DD (double2 x, double2 y, double2 z)
822     {
823     valeur[0] = x;
824     valeur[1] = y;
825     valeur[2] = z;
826     }
827    
828 francois 550 OT_VECTEUR_3DD::OT_VECTEUR_3DD(double2 *xyz1,double2 *xyz2)
829     {
830     valeur[0]=xyz2[0]-xyz1[0];
831     valeur[1]=xyz2[1]-xyz1[1];
832     valeur[2]=xyz2[2]-xyz1[2];
833     }
834    
835 francois 283 OT_VECTEUR_3DD::OT_VECTEUR_3DD (const OT_VECTEUR_3DD & v)
836     {
837     for (int i=0; i<3; i++)
838     valeur[i]=v[i];
839     }
840    
841    
842     double2 & OT_VECTEUR_3DD::operator()(int i)
843     {
844     return valeur[i];
845     }
846    
847     double2 OT_VECTEUR_3DD::operator()(int i) const
848     {
849     return valeur[i];
850     }
851    
852     double2 & OT_VECTEUR_3DD::operator[](int i)
853     {
854     return valeur[i];
855     }
856    
857     double2 OT_VECTEUR_3DD::operator[](int i) const
858     {
859     return valeur[i];
860     }
861    
862     OT_VECTEUR_3DD::operator const double2* () const
863     {
864     return valeur;
865     }
866    
867     OT_VECTEUR_3DD::operator double2* ()
868     {
869     return valeur;
870     }
871    
872     OT_VECTEUR_3DD& OT_VECTEUR_3DD::operator =(const OT_VECTEUR_3DD& v3d)
873     {
874     if (this==&v3d) return *this;
875     valeur [0]=v3d.get_x();
876     valeur [1]=v3d.get_y();
877     valeur [2]=v3d.get_z();
878     return *this;
879     }
880     double2 OT_VECTEUR_3DD::norme()
881     {
882     double2 nor= valeur [0]*valeur [0]+valeur [1]*valeur [1]+valeur [2]*valeur [2];
883     nor=sqrt(nor);
884     return nor;
885     }
886     OT_VECTEUR_3DD OT_VECTEUR_3DD::vecteur_norme()
887     {
888     OT_VECTEUR_3DD nor;
889     double2 zero=0.;
890     double2 norme= this->norme() ;
891     nor=(*this)*(1./norme);
892     return nor;
893     }
894 francois 550 double2 OT_VECTEUR_3DD::get_longueur(void)
895     {
896     return sqrt(valeur[0]*valeur[0]+valeur[1]*valeur[1]+valeur[2]*valeur[2]);
897     }
898 francois 283
899 francois 550 double2 OT_VECTEUR_3DD::get_longueur2(void)
900     {
901     return (valeur[0]*valeur[0]+valeur[1]*valeur[1]+valeur[2]*valeur[2]);
902     }
903    
904 francois 283 double2 OT_VECTEUR_3DD::get_x(void) const
905     {
906     return valeur[0];
907     }
908     double2 & OT_VECTEUR_3DD::x(void)
909     {
910     return valeur[0];
911     }
912     double2 OT_VECTEUR_3DD::get_y(void) const
913     {
914     return valeur[1];
915     }
916     double2 & OT_VECTEUR_3DD::y(void)
917     {
918     return valeur[1];
919     }
920     double2 OT_VECTEUR_3DD::get_z(void) const
921     {
922     return valeur[2];
923     }
924     double2 & OT_VECTEUR_3DD::z(void)
925     {
926     return valeur[2];
927     }
928     void OT_VECTEUR_3DD::change_x(double2 x)
929     {
930     valeur[0] = x;
931     }
932     void OT_VECTEUR_3DD::change_y(double2 y)
933     {
934     valeur[1] = y;
935     }
936     void OT_VECTEUR_3DD::change_z(double2 z)
937     {
938     valeur[2] = z;
939     }
940     double2* OT_VECTEUR_3DD::get_xyz(void)
941     {
942     return valeur;
943     }
944    
945     OT_VECTEUR_3DD OT_VECTEUR_3DD::operator+ (const OT_VECTEUR_3DD & a)
946     {
947     return OT_VECTEUR_3DD(
948     a(0)+valeur[0],a(1)+valeur[1],a(2)+valeur[2]);
949     }
950    
951     OT_VECTEUR_3DD OT_VECTEUR_3DD::operator- (const OT_VECTEUR_3DD & a)
952     {
953     return OT_VECTEUR_3DD(
954     valeur[0]-a(0),valeur[1]-a(1),valeur[2]-a(2));
955     }
956    
957     OT_VECTEUR_3DD operator* (const OT_VECTEUR_3DD & rkV, const double2 a)
958     {
959     return OT_VECTEUR_3DD(
960     a*rkV[0],a*rkV[1],a*rkV[2]);
961     }
962     OT_VECTEUR_3DD operator* (const double2 a, const OT_VECTEUR_3DD & rkV)
963     {
964     return OT_VECTEUR_3DD(
965     a*rkV[0],a*rkV[1],a*rkV[2]);
966     }
967    
968     OT_VECTEUR_3DD operator/ (const OT_VECTEUR_3DD & rkV, const double2 a)
969     {
970     double2 inv = 1/a;
971     return OT_VECTEUR_3DD(
972     inv*rkV[0],inv*rkV[1],inv*rkV[2]);
973     }
974    
975    
976     OT_VECTEUR_3DD operator/ (const double2 a, const OT_VECTEUR_3DD & rkV)
977     {
978     double2 inv = 1/a;
979     return OT_VECTEUR_3DD(
980     inv*rkV[0],inv*rkV[1],inv*rkV[2]);
981     }
982     OT_VECTEUR_3DD OT_VECTEUR_3DD::operator&(const OT_VECTEUR_3DD& v2)
983     {
984     OT_VECTEUR_3DD tmp;
985     tmp[0]= valeur[1]*v2[2]-valeur[2]*v2[1];
986     tmp[1]= valeur[2]*v2[0]-valeur[0]*v2[2];
987     tmp[2]= valeur[0]*v2[1]-valeur[1]*v2[0];
988     return tmp;
989     }
990     double2 OT_VECTEUR_3DD::operator* (const OT_VECTEUR_3DD & a)
991     {
992     return a[0]*valeur[0]+a[1]*valeur[1]+a[2]*valeur[2];
993     }
994    
995    
996     OT_VECTEUR_3DD& OT_VECTEUR_3DD::operator+= (const OT_VECTEUR_3DD& rkV)
997     {
998     valeur[0] = valeur[0] +rkV.valeur[0];
999     valeur[1] = valeur[1] +rkV.valeur[1];
1000     valeur[2] = valeur[2] +rkV.valeur[2];
1001     return *this;
1002     }
1003    
1004     OT_VECTEUR_3DD& OT_VECTEUR_3DD::operator-= (const OT_VECTEUR_3DD& rkV)
1005     {
1006     valeur[0] = valeur[0] -rkV.valeur[0];
1007     valeur[1] = valeur[1] -rkV.valeur[1];
1008     valeur[2] = valeur[2] -rkV.valeur[2];
1009     return *this;
1010     }
1011     OT_VECTEUR_3DD& OT_VECTEUR_3DD::operator*= (double2 fScalar)
1012     {
1013     valeur[0] =valeur[0] * fScalar;
1014     valeur[1] =valeur[1] * fScalar;
1015     valeur[2] =valeur[2] * fScalar;
1016     return *this;
1017     }
1018     OT_VECTEUR_3DD& OT_VECTEUR_3DD::operator/= (double2 fScalar)
1019     {
1020     double2 a = 1/fScalar;
1021     valeur[0] =valeur[0] * a;
1022     valeur[1] =valeur[1] * a;
1023     valeur[2] =valeur[2] * a;
1024     return *this;
1025     }
1026    
1027    
1028     std::ostream & operator << ( std::ostream & __os, const OT_VECTEUR_3DD & __vec)
1029     {
1030     __os << "[ "<<setw(20)<< __vec[0]<<setw(4)<<", "<<setw(20)<<__vec[1]<<", "<<setw(20)<<__vec[2]<<setw(4)<<", "<<setw(20)<<"] ";
1031     return __os;
1032     }
1033    
1034     int OT_VECTEUR_3DD::compare_valeur (const OT_VECTEUR_3DD& mdd) const
1035     {
1036     if (valeur[0]==mdd.valeur[0]&&valeur[1]==mdd.valeur[1] &&valeur[2]==mdd.valeur[2])
1037     return 1;
1038     else return 0;
1039    
1040     }
1041    
1042     bool OT_VECTEUR_3DD::operator==(const OT_VECTEUR_3DD& mdd) const
1043     {
1044     return compare_valeur(mdd);
1045     }
1046    
1047    
1048    
1049    
1050    
1051    
1052    
1053    
1054    
1055    
1056     OT_VECTEUR_4DD::OT_VECTEUR_4DD (double2 v[4])
1057     {
1058     for (int i=0; i<4; i++)
1059     valeur[i]=v[i];
1060     }
1061    
1062     OT_VECTEUR_4DD::OT_VECTEUR_4DD (double2 x, double2 y, double2 z, double2 w)
1063     {
1064     valeur[0] = x;
1065     valeur[1] = y;
1066     valeur[2] = z;
1067     valeur[3] = w;
1068     }
1069    
1070     OT_VECTEUR_4DD::OT_VECTEUR_4DD (const OT_VECTEUR_4DD & v)
1071     {
1072     for (int i=0; i<4; i++)
1073     valeur[i]=v[i];
1074     }
1075    
1076    
1077     double2 & OT_VECTEUR_4DD::operator()(int i)
1078     {
1079     return valeur[i];
1080     }
1081    
1082     double2 OT_VECTEUR_4DD::operator()(int i) const
1083     {
1084     return valeur[i];
1085     }
1086    
1087     double2 & OT_VECTEUR_4DD::operator[](int i)
1088     {
1089     return valeur[i];
1090     }
1091    
1092     double2 OT_VECTEUR_4DD::operator[](int i) const
1093     {
1094     return valeur[i];
1095     }
1096    
1097     OT_VECTEUR_4DD::operator const double2* () const
1098     {
1099     return valeur;
1100     }
1101    
1102     OT_VECTEUR_4DD::operator double2* ()
1103     {
1104     return valeur;
1105     }
1106    
1107     OT_VECTEUR_4DD& OT_VECTEUR_4DD::operator =(const OT_VECTEUR_4DD& v4d)
1108     {
1109     if (this==&v4d) return *this;
1110     valeur [0]=v4d.get_x();
1111     valeur [1]=v4d.get_y();
1112     valeur [2]=v4d.get_z();
1113     valeur [3]=v4d.get_w();
1114     return *this;
1115     }
1116     double2 OT_VECTEUR_4DD::norme()
1117     {
1118     double2 nor= valeur [0]*valeur [0]+valeur [1]*valeur [1]+valeur [2]*valeur [2]+valeur [3]*valeur [3];
1119     nor=nor^0.5;
1120     return nor;
1121     }
1122     OT_VECTEUR_4DD OT_VECTEUR_4DD::vecteur_norme()
1123     {
1124     OT_VECTEUR_4DD nor;
1125     double2 zero=0.;
1126     double2 norme= this->norme() ;
1127     assert(norme!=zero) ;
1128     nor=(*this)*(1./norme);
1129     return nor;
1130     }
1131    
1132 francois 363 OT_VECTEUR_4DD OT_VECTEUR_4DD::vecteur_norme_3d()
1133     {
1134     double2 norme= valeur [0]*valeur [0]+valeur [1]*valeur [1]+valeur [2]*valeur [2];
1135     norme=norme^0.5;
1136     OT_VECTEUR_4DD nor;
1137     double2 zero=0.;
1138     assert(norme!=zero) ;
1139     nor.change_x(valeur[0]/norme);
1140     nor.change_y(valeur[1]/norme);
1141     nor.change_z(valeur[2]/norme);
1142     nor.change_w(valeur[3]);
1143     return nor;
1144     }
1145 francois 283 double2 OT_VECTEUR_4DD::get_x(void) const
1146     {
1147     return valeur[0];
1148     }
1149     double2 & OT_VECTEUR_4DD::x(void)
1150     {
1151     return valeur[0];
1152     }
1153     double2 OT_VECTEUR_4DD::get_y(void) const
1154     {
1155     return valeur[1];
1156     }
1157     double2 & OT_VECTEUR_4DD::y(void)
1158     {
1159     return valeur[1];
1160     }
1161     double2 OT_VECTEUR_4DD::get_z(void) const
1162     {
1163     return valeur[2];
1164     }
1165     double2 & OT_VECTEUR_4DD::z(void)
1166     {
1167     return valeur[2];
1168     }
1169     double2 OT_VECTEUR_4DD::get_w(void) const
1170     {
1171     return valeur[3];
1172     }
1173     double2 & OT_VECTEUR_4DD::w(void)
1174     {
1175     return valeur[3];
1176     }
1177     void OT_VECTEUR_4DD::change_x(double2 x)
1178     {
1179     valeur[0] = x;
1180     }
1181     void OT_VECTEUR_4DD::change_y(double2 y)
1182     {
1183     valeur[1] = y;
1184     }
1185     void OT_VECTEUR_4DD::change_z(double2 z)
1186     {
1187     valeur[2] = z;
1188     }
1189     void OT_VECTEUR_4DD::change_w(double2 w)
1190     {
1191     valeur[3] = w;
1192     }
1193     double2* OT_VECTEUR_4DD::get_xyzw(void)
1194     {
1195     return valeur;
1196     }
1197    
1198     OT_VECTEUR_4DD OT_VECTEUR_4DD::operator+ (const OT_VECTEUR_4DD & a)
1199     {
1200     return OT_VECTEUR_4DD(
1201     a(0)+valeur[0],a(1)+valeur[1],a(2)+valeur[2],a(3)+valeur[3]
1202     );
1203     }
1204    
1205     OT_VECTEUR_4DD OT_VECTEUR_4DD::operator- (const OT_VECTEUR_4DD & a)
1206     {
1207     return OT_VECTEUR_4DD(
1208     valeur[0]-a(0),valeur[1]-a(1),valeur[2]-a(2),valeur[3]-a(3)
1209     );
1210     }
1211    
1212     OT_VECTEUR_4DD operator* (const OT_VECTEUR_4DD & rkV, const double2 a)
1213     {
1214     return OT_VECTEUR_4DD(
1215     a*rkV[0],a*rkV[1],a*rkV[2],a*rkV[3]
1216     );
1217     }
1218     OT_VECTEUR_4DD operator* (const double2 a, const OT_VECTEUR_4DD & rkV)
1219     {
1220     return OT_VECTEUR_4DD(
1221     a*rkV[0],a*rkV[1],a*rkV[2],a*rkV[3]
1222     );
1223     }
1224    
1225     OT_VECTEUR_4DD operator/ (const OT_VECTEUR_4DD & rkV, const double2 a)
1226     {
1227     double2 inv = 1/a;
1228     return OT_VECTEUR_4DD(
1229     inv*rkV[0],inv*rkV[1],inv*rkV[2],inv*rkV[3]
1230     );
1231     }
1232    
1233    
1234     OT_VECTEUR_4DD operator/ (const double2 a, const OT_VECTEUR_4DD & rkV)
1235     {
1236     double2 inv = 1/a;
1237     return OT_VECTEUR_4DD(
1238     inv*rkV[0],inv*rkV[1],inv*rkV[2],inv*rkV[3]
1239     );
1240     }
1241     OT_VECTEUR_4DD OT_VECTEUR_4DD::operator^(const OT_VECTEUR_4DD& v2)
1242     {
1243     OT_VECTEUR_4DD tmp;
1244     tmp[0]= valeur[1]*v2[2]-valeur[2]*v2[1];
1245     tmp[1]= valeur[2]*v2[0]-valeur[0]*v2[2];
1246     tmp[2]= valeur[0]*v2[1]-valeur[1]*v2[0];
1247     tmp[3]= 0;
1248     //tmp[3]= valeur[1]*v2[2]-valeur[2]*v2[1];
1249     return tmp;
1250     }
1251     double2 OT_VECTEUR_4DD::operator* (const OT_VECTEUR_4DD & a)
1252     {
1253     return a[0]*valeur[0]+a[1]*valeur[1]+a[2]*valeur[2]+a[3]*valeur[3];
1254     }
1255    
1256    
1257     OT_VECTEUR_4DD& OT_VECTEUR_4DD::operator+= (const OT_VECTEUR_4DD& rkV)
1258     {
1259     valeur[0] = valeur[0] +rkV.valeur[0];
1260     valeur[1] = valeur[1] +rkV.valeur[1];
1261     valeur[2] = valeur[2] +rkV.valeur[2];
1262     valeur[3] = valeur[3] +rkV.valeur[3];
1263     return *this;
1264     }
1265    
1266     OT_VECTEUR_4DD& OT_VECTEUR_4DD::operator-= (const OT_VECTEUR_4DD& rkV)
1267     {
1268     valeur[0] = valeur[0] -rkV.valeur[0];
1269     valeur[1] = valeur[1] -rkV.valeur[1];
1270     valeur[2] = valeur[2] -rkV.valeur[2];
1271     valeur[3] = valeur[3] -rkV.valeur[3];
1272     return *this;
1273     }
1274     OT_VECTEUR_4DD& OT_VECTEUR_4DD::operator*= (double2 fScalar)
1275     {
1276     valeur[0] =valeur[0] * fScalar;
1277     valeur[1] =valeur[1] * fScalar;
1278     valeur[2] =valeur[2] * fScalar;
1279     valeur[3] =valeur[3] * fScalar;
1280     return *this;
1281     }
1282     OT_VECTEUR_4DD& OT_VECTEUR_4DD::operator/= (double2 fScalar)
1283     {
1284     double2 a = 1/fScalar;
1285     valeur[0] =valeur[0] * a;
1286     valeur[1] =valeur[1] * a;
1287     valeur[2] =valeur[2] * a;
1288     valeur[3] =valeur[3] * a;
1289     return *this;
1290     }
1291    
1292    
1293     std::ostream & operator << ( std::ostream & __os, const OT_VECTEUR_4DD & __vec)
1294     {
1295     __os << "[ "<<setw(20)<< __vec[0]<<setw(4)<<", "<<setw(20)<<__vec[1]<<", "<<setw(20)<<__vec[2]<<setw(4)<<", "<<setw(20)<<__vec[3]<<setw(4)<<"] ";
1296     return __os;
1297     }
1298    
1299     int OT_VECTEUR_4DD::compare_valeur (const OT_VECTEUR_4DD& mdd) const
1300     {
1301     if (valeur[0]==mdd.valeur[0]&&valeur[1]==mdd.valeur[1] &&valeur[2]==mdd.valeur[2]&&valeur[3]==mdd.valeur[3])
1302     return 1;
1303     else return 0;
1304    
1305     }
1306    
1307     bool OT_VECTEUR_4DD::operator==(const OT_VECTEUR_4DD& mdd) const
1308     {
1309     return compare_valeur(mdd);
1310     }
1311    
1312    
1313    
1314 francois 743 bool OT_VECTEUR_4DD::est_nul(void)
1315     {
1316     double2 ZERO(0.);
1317     if (valeur[0]==ZERO&&valeur[1]==ZERO &&valeur[2]==ZERO&&valeur[3]==ZERO)
1318     return true;
1319     return false;
1320     }
1321 francois 283
1322 francois 743 bool OT_VECTEUR_4DD::est_nul_3d(void)
1323     {
1324     double2 ZERO(0.);
1325     if (valeur[0]==ZERO&&valeur[1]==ZERO &&valeur[2]==ZERO)
1326     return true;
1327     return false;
1328     }
1329 francois 283
1330    
1331    
1332    
1333    
1334    
1335    
1336    
1337    
1338    
1339    
1340 francois 743
1341 francois 283 int OPERATEUR::egal(double a,double b,double eps)
1342     {
1343     double res=fabs(a-b);
1344     if (res<eps) return(1);
1345     return(0);
1346     }
1347    
1348 couturad 906 int OPERATEUR::egal(double* xyz1, double* xyz2, double eps)
1349     {
1350     if(fabs(xyz1[0]-xyz2[0])<eps)
1351     if(fabs(xyz1[1]-xyz2[1])<eps)
1352     if(fabs(xyz1[2]-xyz2[2])<eps)
1353     return 1;
1354     return 0;
1355     }
1356 francois 283
1357 couturad 906
1358    
1359 francois 283 double OPERATEUR::qualite_triangle(double* noeud1,double* noeud2,double* noeud3)
1360     {
1361     OT_VECTEUR_3D vec12(noeud1,noeud2);
1362     OT_VECTEUR_3D vec13(noeud1,noeud3);
1363     OT_VECTEUR_3D vec23(noeud2,noeud3);
1364     double d12=vec12.get_longueur();
1365     double d13=vec13.get_longueur();
1366     double d23=vec23.get_longueur();
1367     double dmax=max(d12,d13);
1368     dmax=max(dmax,d23);
1369     double p=0.5*(d12+d13+d23);
1370     double crit=(p-d12)*(p-d13)*(p-d23)/p;
1371     if (crit<0) crit=0.;
1372     crit=sqrt(crit);
1373     crit=crit/dmax;
1374     crit=crit*3.4641101614;
1375     return crit;
1376     }
1377    
1378     void OPERATEUR::doubleto2int(double val,int& val1,int& val2)
1379     {
1380     int *p=(int*)(&(val));
1381     val1=(*p);
1382     val2=(*(p+1));
1383     }
1384    
1385     double OPERATEUR::qualite_tetra(double* noeud1,double* noeud2,double* noeud3,double *noeud4)
1386     {
1387     OT_VECTEUR_3D ab(noeud1,noeud2);
1388     OT_VECTEUR_3D ac(noeud1,noeud3);
1389     OT_VECTEUR_3D cb(noeud3,noeud2);
1390     OT_VECTEUR_3D db(noeud4,noeud2);
1391     OT_VECTEUR_3D da(noeud4,noeud1);
1392     OT_VECTEUR_3D dc(noeud4,noeud3);
1393     double dab=ab.get_longueur();
1394     double dac=ac.get_longueur();
1395     double dcb=cb.get_longueur();
1396     double ddb=db.get_longueur();
1397     double dda=da.get_longueur();
1398     double ddc=dc.get_longueur();
1399    
1400     double dmax=max(dab,dac);
1401     dmax=max(dmax,dcb);
1402     dmax=max(dmax,ddb);
1403     dmax=max(dmax,dda);
1404     dmax=max(dmax,ddc);
1405    
1406     OT_VECTEUR_3D pvec=ab&ac;
1407     double vol=pvec*da;
1408     if (vol>0.) return 0.;
1409     vol=fabs(vol);
1410     double som=pvec.get_longueur()/2.;
1411     pvec=ab&da;
1412     som=som+pvec.get_longueur()/2.;
1413     pvec=ac&dc;
1414     som=som+pvec.get_longueur()/2.;
1415     pvec=cb&db;
1416     som=som+pvec.get_longueur()/2.;
1417     double crit=vol/som/dmax/0.408249;
1418     return crit;
1419     }
1420    
1421 francois 418 int OPERATEUR::estdansletetra(double *xyz1,double *xyz2,double *xyz3,double *xyz4,double x,double y, double z)
1422     {
1423     OT_VECTEUR_3D v1(xyz2[0]-xyz1[0],xyz2[1]-xyz1[1],xyz2[2]-xyz1[2]);
1424     OT_VECTEUR_3D v2(xyz3[0]-xyz1[0],xyz3[1]-xyz1[1],xyz3[2]-xyz1[2]);
1425     OT_VECTEUR_3D v3(xyz4[0]-xyz1[0],xyz4[1]-xyz1[1],xyz4[2]-xyz1[2]);
1426     OT_VECTEUR_3D v4(x-xyz1[0],y-xyz1[1],z-xyz1[2]);
1427     OT_MATRICE_3D mat(v1,v2,v3);
1428     OT_MATRICE_3D mat1(v4,v2,v3);
1429     OT_MATRICE_3D mat2(v1,v4,v3);
1430     OT_MATRICE_3D mat3(v1,v2,v4);
1431     double det=mat.get_determinant();
1432     double xsi=mat1.get_determinant()/det;
1433     double eta=mat2.get_determinant()/det;
1434     double dseta=mat3.get_determinant()/det;
1435     int reponse1=1;
1436     int reponse2=1;
1437     double eps=0.0000001;
1438     if (xsi<-eps) reponse1=0;
1439     if (eta<-eps) reponse1=0;
1440     if (dseta<-eps) reponse1=0;
1441     if (xsi+eta+dseta>1.+eps) reponse1=0;
1442     if (xsi<eps) reponse2=0;
1443     if (eta<eps) reponse2=0;
1444     if (dseta<eps) reponse2=0;
1445     if (xsi+eta+dseta>1.-eps) reponse2=0;
1446     int face1=0,face2=0,face3=0,face4=0;
1447     if ((reponse1==1) && (reponse2==0))
1448     {
1449     if (dseta<eps) face1=1;
1450     if (eta<eps) face2=1;
1451     if (xsi<eps) face3=1;
1452     if (xsi+eta+dseta>1.-eps) face4=1;
1453     }
1454     return reponse1+2*reponse2+4*face1+8*face2+16*face3+32*face4;
1455     }
1456 francois 283
1457    
1458 francois 418
1459    
1460     int OPERATEUR::projeteestdansletriangle(double *xyz1,double *xyz2,double *xyz3,double x,double y, double z)
1461     {
1462     OT_VECTEUR_3D vec1(xyz1,xyz2);
1463     OT_VECTEUR_3D vec2(xyz1,xyz3);
1464     OT_VECTEUR_3D n=vec1&vec2;
1465     n.norme();
1466     double d=-(n.get_x()*xyz1[0]+n.get_y()*xyz1[1]+n.get_z()*xyz1[2]);
1467     double t=-(n.get_x()*x+n.get_y()*y+n.get_z()*z+d);
1468     double xx=x+t*n.get_x();
1469     double yy=y+t*n.get_y();
1470     double zz=z+t*n.get_z();
1471     return estdansletriangle(xyz1,xyz2,xyz3,xx,yy,zz);
1472     }
1473    
1474    
1475     int OPERATEUR::estdansletriangle(double *xyz1,double *xyz2,double *xyz3,double x,double y, double z)
1476     {
1477     OT_VECTEUR_3D vec1(xyz1,xyz2);
1478     OT_VECTEUR_3D vec2(xyz1,xyz3);
1479     OT_VECTEUR_3D vecb(x-xyz1[0],y-xyz1[1],z-xyz1[2]);
1480     OT_MATRICE_3D m(vec1,vec2,vecb);
1481     double det=m.get_determinant();
1482     double eps=1e-10*((vec1&vec2)*(vec1&vec2));
1483     if (fabs(det)>eps) return 0;
1484     double ab[3]={vec1.get_x(),vec1.get_y(),vec1.get_z()};
1485     double ac[3]={vec2.get_x(),vec2.get_y(),vec2.get_z()};
1486     double b[3]={vecb.get_x(),vecb.get_y(),vecb.get_z()};
1487     det=ab[0]*ac[1]-ab[1]*ac[0];
1488     int n1,n2;
1489     if (fabs(det)>eps)
1490     {
1491     n1=0;
1492     n2=1;
1493     }
1494     else
1495     {
1496     det=ab[0]*ac[2]-ab[2]*ac[0];
1497     if (fabs(det)>eps)
1498     {
1499     n1=0;
1500     n2=2;
1501     }
1502     else
1503     {
1504     det=ab[1]*ac[2]-ab[2]*ac[1];
1505     n1=1;
1506     n2=2;
1507     }
1508     }
1509     double xsi=(b[n1]*ac[n2]-b[n2]*ac[n1])/det;
1510     double eta=(ab[n1]*b[n2]-ab[n2]*b[n1])/det;
1511     int reponse1=1;
1512     int reponse2=1;
1513     eps=0.0000001;
1514     if (xsi<-eps) reponse1=0;
1515     if (eta<-eps) reponse1=0;
1516     if (xsi+eta>1.+eps) reponse1=0;
1517     if (xsi<eps) reponse2=0;
1518     if (eta<eps) reponse2=0;
1519     if (xsi+eta>1.-eps) reponse2=0;
1520     int arete1=0,arete2=0,arete3=0;
1521     if ((reponse1==1) && (reponse2==0))
1522     {
1523     if (eta<eps) arete1=1;
1524     if (xsi<eps) arete3=1;
1525     if (xsi+eta>1.-eps) arete2=1;
1526     }
1527     return reponse1+2*reponse2+4*arete1+8*arete2+16*arete3;
1528     }
1529    
1530    
1531    
1532    
1533     int OPERATEUR::compare_etat_tetra(int etat,int valeur)
1534     {
1535     int faceinte=etat&1;
1536     int inte=(etat&2)>>1;
1537     int face1=(etat&4)>>2;
1538     int face2=(etat&8)>>3;
1539     int face3=(etat&16)>>4;
1540     int face4=(etat&32)>>5;
1541     if ((faceinte==1) && (valeur==INTERIEUR)) return 1;
1542     if ((inte==1) && (valeur==INTERIEUR)) return 1;
1543     if (/*(face==0) && */(inte==1) && (valeur==STRICTINTERIEUR)) return 1;
1544     if ((faceinte==1) && (inte==0) && (valeur==SUR_FACE)) return 1;
1545     if ((face1==1) && (valeur==FACE1)) return 1;
1546     if ((face2==1) && (valeur==FACE2)) return 1;
1547     if ((face3==1) && (valeur==FACE3)) return 1;
1548     if ((face4==1) && (valeur==FACE4)) return 1;
1549     if ((face1==1) && (face2==1) && (valeur==ARETE1)) return 1;
1550     if ((face1==1) && (face3==1) && (valeur==ARETE2)) return 1;
1551     if ((face2==1) && (face3==1) && (valeur==ARETE3)) return 1;
1552     if ((face1==1) && (face4==1) && (valeur==ARETE4)) return 1;
1553     if ((face2==1) && (face4==1) && (valeur==ARETE5)) return 1;
1554     if ((face3==1) && (face4==1) && (valeur==ARETE6)) return 1;
1555     if ((face1==1) && (face2==1) && (face3==1) && (valeur==SOMMET1)) return 1;
1556     if ((face1==1) && (face2==1) && (face4==1) && (valeur==SOMMET2)) return 1;
1557     if ((face1==1) && (face3==1) && (face4==1) && (valeur==SOMMET3)) return 1;
1558     if ((face2==1) && (face3==1) && (face4==1) && (valeur==SOMMET4)) return 1;
1559     return 0;
1560     }
1561    
1562     int OPERATEUR::compare_etat_triangle(int etat,int valeur)
1563     {
1564     int areteinte=etat&1;
1565     int inte=(etat&2)>>1;
1566     int arete1=(etat&4)>>2;
1567     int arete2=(etat&8)>>3;
1568     int arete3=(etat&16)>>4;
1569     if ((areteinte==1) && (valeur==INTERIEUR)) return 1;
1570     if ((inte==1) && (valeur==INTERIEUR)) return 1;
1571     if ((inte==1) && (valeur==STRICTINTERIEUR)) return 1;
1572     if ((areteinte==1) && (inte==0) && (valeur==SUR_ARETE)) return 1;
1573     if ((arete1==1) && (valeur==ARETE1)) return 1;
1574     if ((arete2==1) && (valeur==ARETE2)) return 1;
1575     if ((arete3==1) && (valeur==ARETE3)) return 1;
1576     if ((arete1==1) && (arete2==1) && (valeur==SOMMET2)) return 1;
1577     if ((arete2==1) && (arete3==1) && (valeur==SOMMET3)) return 1;
1578     if ((arete3==1) && (arete1==1) && (valeur==SOMMET1)) return 1;
1579     return 0;
1580     }
1581    
1582    
1583 francois 283 template <int N> DOUBLEN<N>::DOUBLEN()
1584     {
1585     tab[0]=-1;
1586     }
1587    
1588     template <int N> DOUBLEN<N>::DOUBLEN(DOUBLEN& mdd)
1589     {
1590     for (int i=0;i<N;i++)
1591     tab[i]=mdd.tab[i];
1592     }
1593    
1594     template <int N> DOUBLEN<N>::DOUBLEN(double *v)
1595     {
1596     for (int i=0;i<N;i++)
1597     tab[i]=v[i];
1598     }
1599    
1600     template <int N> DOUBLEN<N>::~DOUBLEN()
1601     {
1602     }
1603    
1604     template <int N> double DOUBLEN<N>::get_valeur(int num)
1605     {
1606     return tab[num];
1607     }
1608    
1609     template <int N> void DOUBLEN<N>::change_valeur(int num,double val)
1610     {
1611     tab[num]=val;
1612     }
1613    
1614     template <int N> DOUBLEN<N> & DOUBLEN<N>::operator=(DOUBLEN<N> & mdd)
1615     {
1616     for (int i=0;i<N;i++)
1617     tab[i]=mdd.tab[i];
1618     return *this;
1619     }
1620    
1621     template <int N> DOUBLEN<N> & DOUBLEN<N>::operator=(const DOUBLEN<N> & mdd)
1622     {
1623     for (int i=0;i<N;i++)
1624     tab[i]=mdd.tab[i];
1625     return *this;
1626     }
1627    
1628    
1629     #ifdef BORLANDCPP
1630     // instanciation pour Windows (Borland c++)
1631     template class __export DOUBLEN<1>;
1632     template class __export DOUBLEN<4>;
1633     template class __export DOUBLEN<10>;
1634     #else
1635     #ifdef GCC
1636     // instanciation pour Unix (g++)
1637     template class DOUBLEN<1>;
1638     template class DOUBLEN<4>;
1639     template class DOUBLEN<10>;
1640     #else
1641     // instanciation pour Windows (Visual Studio)
1642     #endif
1643     #endif
1644    
1645