ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/outil/src/ot_mathematique.cpp
Revision: 743
Committed: Mon Oct 5 15:00:32 2015 UTC (9 years, 7 months ago) by francois
File size: 35406 byte(s)
Log Message:
Suppression des vecteurs nuls dans la vectorisation ca devient applicable pour une sphere complete

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