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, 9 months ago) by couturad
File size: 36654 byte(s)
Log Message:
Nouveau opencascade commit 1

File Contents

# Content
1 //------------------------------------------------------------
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 OT_QUATERNION::OT_QUATERNION(void)
36 {
37 valeur[0]=0;
38 valeur[1]=0;
39 valeur[2]=0;
40 valeur[3]=0;
41 }
42
43 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 OT_VECTEUR_3D::OT_VECTEUR_3D(double x,double y,double z)
110 {
111 valeur[0]=x;
112 valeur[1]=y;
113 valeur[2]=z;
114 }
115
116
117 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
124 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
131 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 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 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 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
899 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 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 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 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 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
1322 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
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341 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 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
1357
1358
1359 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 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
1457
1458
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 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