ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/dly_tetra.cpp
Revision: 1158
Committed: Thu Jun 13 22:18:49 2024 UTC (10 months, 4 weeks ago) by francois
File size: 15543 byte(s)
Log Message:
compatibilité Ubuntu 22.04
Suppression des refeences à Windows
Ajout d'une banière

File Contents

# User Rev Content
1 francois 1158 //####//------------------------------------------------------------
2     //####//------------------------------------------------------------
3     //####// MAGiC
4     //####// Jean Christophe Cuilliere et Vincent FRANCOIS
5     //####// Departement de Genie Mecanique - UQTR
6     //####//------------------------------------------------------------
7     //####// MAGIC est un projet de recherche de l equipe ERICCA
8     //####// du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
9     //####// http://www.uqtr.ca/ericca
10     //####// http://www.uqtr.ca/
11     //####//------------------------------------------------------------
12     //####//------------------------------------------------------------
13     //####//
14     //####// dly_tetra.cpp
15     //####//
16     //####//------------------------------------------------------------
17     //####//------------------------------------------------------------
18     //####// COPYRIGHT 2000-2024
19     //####// jeu 13 jun 2024 11:58:55 EDT
20     //####//------------------------------------------------------------
21     //####//------------------------------------------------------------
22 francois 285
23    
24     #include "gestionversion.h"
25    
26 francois 313 #include <stdlib.h>
27 francois 285 #include "dly_noeud.h"
28     #include "dly_segment.h"
29     #include "dly_tetra.h"
30     #include "dly_triangle.h"
31     #include "ot_mathematique.h"
32     #include "fct_taille.h"
33 francois 1158 #include "robustpredicates.hxx"
34 francois 285 #include <math.h>
35     #include <string.h>
36    
37    
38    
39     DLY_TETRA::DLY_TETRA(class DLY_NOEUD *no1,class DLY_NOEUD *no2,class DLY_NOEUD *no3,class DLY_NOEUD *no4):noeud1(no1),noeud2(no2),noeud3(no3),noeud4(no4),feuille(1),voisin1(NULL),voisin2(NULL),voisin3(NULL),voisin4(NULL),normal1(NULL),normal2(NULL),normal3(NULL),normal4(NULL),qualite(-100.)
40     {
41     }
42    
43     DLY_TETRA::DLY_TETRA(class DLY_NOEUD *no1,class DLY_NOEUD *no2,class DLY_NOEUD *no3,class DLY_NOEUD *no4,double qual):noeud1(no1),noeud2(no2),noeud3(no3),noeud4(no4),feuille(1),voisin1(NULL),voisin2(NULL),voisin3(NULL),voisin4(NULL),normal1(NULL),normal2(NULL),normal3(NULL),normal4(NULL),qualite(qual)
44     {
45     }
46     DLY_TETRA::DLY_TETRA(DLY_TETRA& mdd):noeud1(mdd.noeud1),noeud2(mdd.noeud2),noeud3(mdd.noeud3),noeud4(mdd.noeud4),feuille(mdd.feuille),voisin1(mdd.voisin1),voisin2(mdd.voisin2),voisin3(mdd.voisin3),voisin4(mdd.voisin4),qualite(mdd.qualite)
47     {
48     if (mdd.normal1==NULL) normal1=NULL;
49     else normal1=new double[3];
50     if (mdd.normal2==NULL) normal2=NULL;
51     else normal2=new double[3];
52     if (mdd.normal3==NULL) normal3=NULL;
53     else normal3=new double[3];
54     if (mdd.normal4==NULL) normal4=NULL;
55     else normal4=new double[3];
56     if (mdd.normal1!=NULL) memcpy(normal1,mdd.normal1,3*sizeof(double));
57     if (mdd.normal2!=NULL) memcpy(normal2,mdd.normal2,3*sizeof(double));
58     if (mdd.normal3!=NULL) memcpy(normal3,mdd.normal3,3*sizeof(double));
59     if (mdd.normal4!=NULL) memcpy(normal4,mdd.normal4,3*sizeof(double));
60     }
61     DLY_TETRA::~DLY_TETRA()
62     {
63     if (normal1!=NULL) delete [] normal1;
64     if (normal2!=NULL) delete [] normal2;
65     if (normal3!=NULL) delete [] normal3;
66     if (normal4!=NULL) delete [] normal4;
67     }
68    
69     DLY_NOEUD* DLY_TETRA::get_noeud1(void)
70     {
71     return noeud1;
72     }
73     DLY_NOEUD* DLY_TETRA::get_noeud2(void)
74     {
75     return noeud2;
76     }
77     DLY_NOEUD* DLY_TETRA::get_noeud3(void)
78     {
79     return noeud3;
80     }
81     DLY_NOEUD* DLY_TETRA::get_noeud4(void)
82     {
83     return noeud4;
84     }
85    
86     int DLY_TETRA::get_feuille(void)
87     {
88     return feuille;
89     }
90    
91     void DLY_TETRA::change_feuille(int num)
92     {
93     feuille=num;
94     }
95    
96     DLY_TRIANGLE* DLY_TETRA::get_triangle1(void)
97     {
98     return tri1;
99     }
100     DLY_TRIANGLE* DLY_TETRA::get_triangle2(void)
101     {
102     return tri2;
103     }
104     DLY_TRIANGLE* DLY_TETRA::get_triangle3(void)
105     {
106     return tri3;
107     }
108     DLY_TRIANGLE* DLY_TETRA::get_triangle4(void)
109     {
110     return tri4;
111     }
112     DLY_SEGMENT* DLY_TETRA::get_segment1(void)
113     {
114     return seg1;
115     }
116     DLY_SEGMENT* DLY_TETRA::get_segment2(void)
117     {
118     return seg2;
119     }
120     DLY_SEGMENT* DLY_TETRA::get_segment3(void)
121     {
122     return seg3;
123     }
124     DLY_SEGMENT* DLY_TETRA::get_segment4(void)
125     {
126     return seg4;
127     }
128     DLY_SEGMENT* DLY_TETRA::get_segment5(void)
129     {
130     return seg5;
131     }
132     DLY_SEGMENT* DLY_TETRA::get_segment6(void)
133     {
134     return seg6;
135     }
136     void DLY_TETRA::change_triangle1(DLY_TRIANGLE* tri)
137     {
138     tri1=tri;
139     tri1->insere_lien_tetra(this);
140     }
141     void DLY_TETRA::change_triangle2(DLY_TRIANGLE* tri)
142     {
143     tri2=tri;
144     tri2->insere_lien_tetra(this);
145     }
146     void DLY_TETRA::change_triangle3(DLY_TRIANGLE* tri)
147     {
148     tri3=tri;
149     tri3->insere_lien_tetra(this);
150     }
151     void DLY_TETRA::change_triangle4(DLY_TRIANGLE* tri)
152     {
153     tri4=tri;
154     tri4->insere_lien_tetra(this);
155     }
156     void DLY_TETRA::change_segment1(DLY_SEGMENT* seg)
157     {
158     seg1=seg;
159     seg1->insere_lien_tetra(this);
160     }
161     void DLY_TETRA::change_segment2(DLY_SEGMENT* seg)
162     {
163     seg2=seg;
164     seg2->insere_lien_tetra(this);
165     }
166     void DLY_TETRA::change_segment3(DLY_SEGMENT* seg)
167     {
168     seg3=seg;
169     seg3->insere_lien_tetra(this);
170     }
171     void DLY_TETRA::change_segment4(DLY_SEGMENT* seg)
172     {
173     seg4=seg;
174     seg4->insere_lien_tetra(this);
175     }
176     void DLY_TETRA::change_segment5(DLY_SEGMENT* seg)
177     {
178     seg5=seg;
179     seg5->insere_lien_tetra(this);
180     }
181     void DLY_TETRA::change_segment6(DLY_SEGMENT* seg)
182     {
183     seg6=seg;
184     seg6->insere_lien_tetra(this);
185     }
186    
187     DLY_TETRA* DLY_TETRA::get_voisin1(void)
188     {
189     return voisin1;
190     }
191     DLY_TETRA* DLY_TETRA::get_voisin2(void)
192     {
193     return voisin2;
194     }
195     DLY_TETRA* DLY_TETRA::get_voisin3(void)
196     {
197     return voisin3;
198     }
199     DLY_TETRA* DLY_TETRA::get_voisin4(void)
200     {
201     return voisin4;
202     }
203     void DLY_TETRA::change_voisin1(DLY_TETRA* tet)
204     {
205     voisin1=tet;
206     }
207     void DLY_TETRA::change_voisin2(DLY_TETRA* tet)
208     {
209     voisin2=tet;
210     }
211     void DLY_TETRA::change_voisin3(DLY_TETRA* tet)
212     {
213     voisin3=tet;
214     }
215     void DLY_TETRA::change_voisin4(DLY_TETRA* tet)
216     {
217     voisin4=tet;
218     }
219     void DLY_TETRA::ajoute_fils(DLY_TETRA* tet)
220     {
221     feuille=0;
222     fils.push_back(tet);
223     }
224     void DLY_TETRA::ajoute_fils(std::vector<DLY_TETRA*> &liste)
225     {
226     if (liste.size()==0) return;
227     feuille=0;
228     for (int i=0;i<liste.size();i++)
229     fils.push_back(liste[i]);
230     }
231     int DLY_TETRA::get_nb_fils(void)
232     {
233     return fils.size();
234     }
235     DLY_TETRA* DLY_TETRA::get_fils(int i)
236     {
237     return fils[i];
238     }
239    
240     int DLY_TETRA::point_dans_la_sphere(double x,double y,double z)
241     {
242     double xyz1[3],xyz2[3],xyz3[3],xyz4[3];
243     noeud1->get_coord(xyz1);
244     noeud2->get_coord(xyz2);
245     noeud3->get_coord(xyz3);
246     noeud4->get_coord(xyz4);
247     double xyz[3]={x,y,z};
248     double val=robustPredicates::insphere(xyz1,xyz2,xyz3,xyz4,xyz)*robustPredicates::orient3d(xyz1,xyz2,xyz3,xyz4);
249     return (val > 0) ? 1 : 0;
250     }
251     int DLY_TETRA::point_dans_le_tetra(double x,double y,double z)
252     {
253     double xyz1[3],xyz2[3],xyz3[3],xyz4[3];
254     noeud1->get_coord(xyz1);
255     noeud2->get_coord(xyz2);
256     noeud3->get_coord(xyz3);
257     noeud4->get_coord(xyz4);
258     OT_VECTEUR_3D v1(xyz2[0]-xyz1[0],xyz2[1]-xyz1[1],xyz2[2]-xyz1[2]);
259     OT_VECTEUR_3D v2(xyz3[0]-xyz1[0],xyz3[1]-xyz1[1],xyz3[2]-xyz1[2]);
260     OT_VECTEUR_3D v3(xyz4[0]-xyz1[0],xyz4[1]-xyz1[1],xyz4[2]-xyz1[2]);
261     OT_VECTEUR_3D v4(x-xyz1[0],y-xyz1[1],z-xyz1[2]);
262     OT_MATRICE_3D mat(v1,v2,v3);
263     OT_MATRICE_3D mat1(v4,v2,v3);
264     OT_MATRICE_3D mat2(v1,v4,v3);
265     OT_MATRICE_3D mat3(v1,v2,v4);
266     double det=mat.get_determinant();
267     double xsi=mat1.get_determinant()/det;
268     double eta=mat2.get_determinant()/det;
269     double dseta=mat3.get_determinant()/det;
270     int reponse1=1;
271     double eps=0.0000001;
272     if (xsi<-eps) reponse1=0;
273     if (eta<-eps) reponse1=0;
274     if (dseta<-eps) reponse1=0;
275     if (xsi+eta+dseta>1.+eps) reponse1=0;
276     return reponse1;
277     }
278     double DLY_TETRA::get_volume(void)
279     {
280     return volume;
281     }
282     double DLY_TETRA::calcul_volume(FCT_TAILLE *metrique,int conserve)
283     {
284     double vol;
285     double xyz1[3],xyz2[3],xyz3[3],xyz4[3];
286     noeud1->get_coord(xyz1);
287     noeud2->get_coord(xyz2);
288     noeud3->get_coord(xyz3);
289     noeud4->get_coord(xyz4);
290     OT_VECTEUR_3D vec1(xyz1,xyz2);
291     OT_VECTEUR_3D vec2(xyz1,xyz3);
292     OT_VECTEUR_3D vec3(xyz1,xyz4);
293     double facteur=(vec1&vec2)*vec3;
294     if (metrique==NULL)
295     {
296     vol=facteur/6.;
297     if (conserve) volume=vol;
298     return vol;
299     }
300     static int r=15;
301     static double gaussx[15],gaussy[15],gaussz[15],wi[15];
302     static int first=0;
303     if (first==0)
304     {
305     first=1;
306     double a=0.25;
307     double b2=0.319793627;
308     double b1=0.091971078;
309     double c1=0.724086765;
310     double c2=0.040619116;
311     double d=0.056350832;
312     double e=0.443649167;
313     double w0=0.019753086;
314     double w1=0.011989513;
315     double w2=0.011511367;
316     double w3=0.008818342;
317     gaussx[0]=a;
318     gaussy[0]=a;
319     gaussz[0]=a;
320     wi[0]=w0;
321     gaussx[1]=b1;
322     gaussy[1]=b1;
323     gaussz[1]=b1;
324     wi[1]=w1;
325     gaussx[2]=b1;
326     gaussy[2]=b1;
327     gaussz[2]=c1;
328     wi[2]=w1;
329     gaussx[3]=b1;
330     gaussy[3]=c1;
331     gaussz[3]=b1;
332     wi[3]=w1;
333     gaussx[4]=c1;
334     gaussy[4]=b1;
335     gaussz[4]=b1;
336     wi[4]=w1;
337     gaussx[5]=b2;
338     gaussy[5]=b2;
339     gaussz[5]=b2;
340     wi[5]=w2;
341     gaussx[6]=b2;
342     gaussy[6]=b2;
343     gaussz[6]=c2;
344     wi[6]=w2;
345     gaussx[7]=b2;
346     gaussy[7]=c2;
347     gaussz[7]=b2;
348     wi[7]=w2;
349     gaussx[8]=c2;
350     gaussy[8]=b2;
351     gaussz[8]=b2;
352     wi[8]=w2;
353     gaussx[9]=d;
354     gaussy[9]=d;
355     gaussz[9]=e;
356     wi[9]=w3;
357     gaussx[10]=d;
358     gaussy[10]=e;
359     gaussz[10]=d;
360     wi[10]=w3;
361     gaussx[11]=e;
362     gaussy[11]=d;
363     gaussz[11]=d;
364     wi[11]=w3;
365     gaussx[12]=d;
366     gaussy[12]=e;
367     gaussz[12]=e;
368     wi[12]=w3;
369     gaussx[13]=e;
370     gaussy[13]=d;
371     gaussz[13]=e;
372     wi[13]=w3;
373     gaussx[14]=e;
374     gaussy[14]=e;
375     gaussz[14]=d;
376     wi[14]=w3;
377     }
378     vol=0.;
379     for (int i=0;i<r;i++)
380     {
381     double xyz[3];
382     xyz[0]=(1-gaussx[i]-gaussy[i]-gaussz[i])*xyz1[0]+gaussx[i]*xyz2[0]+gaussy[i]*xyz3[0]+gaussz[i]*xyz4[0];
383     xyz[1]=(1-gaussx[i]-gaussy[i]-gaussz[i])*xyz1[1]+gaussx[i]*xyz2[1]+gaussy[i]*xyz3[1]+gaussz[i]*xyz4[1];
384     xyz[2]=(1-gaussx[i]-gaussy[i]-gaussz[i])*xyz1[2]+gaussx[i]*xyz2[2]+gaussy[i]*xyz3[2]+gaussz[i]*xyz4[2];
385     double tenseur[9];
386     metrique->evaluer(xyz,tenseur);
387     OT_VECTEUR_3D vec1(tenseur[0],tenseur[3],tenseur[6]);
388     OT_VECTEUR_3D vec2(tenseur[1],tenseur[4],tenseur[7]);
389     OT_VECTEUR_3D vec3(tenseur[2],tenseur[5],tenseur[8]);
390     OT_MATRICE_3D mat(vec1,vec2,vec3);
391     vol=vol+sqrt(mat.get_determinant())*wi[i]*facteur;
392     }
393     if (conserve) volume=vol;
394     return vol;
395     }
396    
397     void DLY_TETRA::get_normal1(double *nor)
398     {
399     if (normal1==NULL)
400     {
401     normal1=new double[3];
402     double xyz1[3],xyz2[3],xyz3[3];
403     noeud1->get_coord(xyz1);
404     noeud2->get_coord(xyz2);
405     noeud3->get_coord(xyz3);
406     OT_VECTEUR_3D vec1(xyz1,xyz3);
407     OT_VECTEUR_3D vec2(xyz1,xyz2);
408     OT_VECTEUR_3D n=vec1&vec2;
409     n.norme();
410     normal1[0]=n.get_x();
411     normal1[1]=n.get_y();
412     normal1[2]=n.get_z();
413     }
414     nor[0]=normal1[0];
415     nor[1]=normal1[1];
416     nor[2]=normal1[2];
417     }
418     void DLY_TETRA::get_normal2(double *nor)
419     {
420     if (normal2==NULL)
421     {
422     normal2=new double[3];
423     double xyz1[3],xyz2[3],xyz4[3];
424     noeud1->get_coord(xyz1);
425     noeud2->get_coord(xyz2);
426     noeud4->get_coord(xyz4);
427     OT_VECTEUR_3D vec1(xyz1,xyz2);
428     OT_VECTEUR_3D vec2(xyz1,xyz4);
429     OT_VECTEUR_3D n=vec1&vec2;
430     n.norme();
431     normal2[0]=n.get_x();
432     normal2[1]=n.get_y();
433     normal2[2]=n.get_z();
434     }
435     nor[0]=normal2[0];
436     nor[1]=normal2[1];
437     nor[2]=normal2[2];
438     }
439     void DLY_TETRA::get_normal3(double *nor)
440     {
441     if (normal3==NULL)
442     {
443     normal3=new double[3];
444     double xyz1[3],xyz3[3],xyz4[3];
445     noeud1->get_coord(xyz1);
446     noeud3->get_coord(xyz3);
447     noeud4->get_coord(xyz4);
448     OT_VECTEUR_3D vec1(xyz1,xyz4);
449     OT_VECTEUR_3D vec2(xyz1,xyz3);
450     OT_VECTEUR_3D n=vec1&vec2;
451     n.norme();
452     normal3[0]=n.get_x();
453     normal3[1]=n.get_y();
454     normal3[2]=n.get_z();
455     }
456     nor[0]=normal3[0];
457     nor[1]=normal3[1];
458     nor[2]=normal3[2];
459     }
460     void DLY_TETRA::get_normal4(double *nor)
461     {
462     if (normal4==NULL)
463     {
464     normal4=new double[4];
465     double xyz2[3],xyz3[3],xyz4[3];
466     noeud2->get_coord(xyz2);
467     noeud3->get_coord(xyz3);
468     noeud4->get_coord(xyz4);
469     OT_VECTEUR_3D vec1(xyz2,xyz3);
470     OT_VECTEUR_3D vec2(xyz2,xyz4);
471     OT_VECTEUR_3D n=vec1&vec2;
472     n.norme();
473     normal4[0]=n.get_x();
474     normal4[1]=n.get_y();
475     normal4[2]=n.get_z();
476     }
477     nor[0]=normal4[0];
478     nor[1]=normal4[1];
479     nor[2]=normal4[2];
480     }
481     double DLY_TETRA::get_qualite(void)
482     {
483     if (qualite<-1.)
484     {
485     double xyz1[3],xyz2[3],xyz3[3],xyz4[3];
486     noeud1->get_coord(xyz1);
487     noeud2->get_coord(xyz2);
488     noeud3->get_coord(xyz3);
489     noeud4->get_coord(xyz4);
490     qualite=OPERATEUR::qualite_tetra(xyz1,xyz2,xyz3,xyz4);
491     }
492     return qualite;
493     }
494    
495     double DLY_TETRA::get_critere_subdivision(void)
496     {
497     return get_volume();
498     }
499    
500     void DLY_TETRA::decoupe_noeud(std::vector<double> &lstpt) // a adapter pour taille variable
501     {
502     double xyz1[3],xyz2[3],xyz3[3],xyz4[3],x,y,z;
503     noeud1->get_coord(xyz1);
504     noeud2->get_coord(xyz2);
505     noeud3->get_coord(xyz3);
506     noeud4->get_coord(xyz4);
507     double maxlon;
508     OT_VECTEUR_3D vec1(xyz1,xyz2);
509     double l1=vec1.get_longueur();
510     maxlon=l1;
511     OT_VECTEUR_3D vec2(xyz1,xyz3);
512     double l2=vec2.get_longueur();
513     if (l2>maxlon) maxlon=l2;
514     OT_VECTEUR_3D vec3(xyz1,xyz4);
515     double l3=vec1.get_longueur();
516     if (l3>maxlon) maxlon=l3;
517     OT_VECTEUR_3D vec4(xyz2,xyz3);
518     double l4=vec1.get_longueur();
519     if (l4>maxlon) maxlon=l4;
520     OT_VECTEUR_3D vec5(xyz2,xyz4);
521     double l5=vec1.get_longueur();
522     if (l5>maxlon) maxlon=l5;
523     OT_VECTEUR_3D vec6(xyz3,xyz4);
524     double l6=vec1.get_longueur();
525     if (l6>maxlon) maxlon=l6;
526     {
527     x=0.5*(xyz1[0]+xyz2[0]);
528     y=0.5*(xyz1[1]+xyz2[1]);
529     z=0.5*(xyz1[2]+xyz2[2]);
530     }
531     lstpt.push_back(x);
532     lstpt.push_back(y);
533     lstpt.push_back(z);
534     {
535     x=0.5*(xyz1[0]+xyz3[0]);
536     y=0.5*(xyz1[1]+xyz3[1]);
537     z=0.5*(xyz1[2]+xyz3[2]);
538     }
539     lstpt.push_back(x);
540     lstpt.push_back(y);
541     lstpt.push_back(z);
542     {
543     x=0.5*(xyz1[0]+xyz4[0]);
544     y=0.5*(xyz1[1]+xyz4[1]);
545     z=0.5*(xyz1[2]+xyz4[2]);
546     }
547     lstpt.push_back(x);
548     lstpt.push_back(y);
549     lstpt.push_back(z);
550     {
551     x=0.5*(xyz2[0]+xyz3[0]);
552     y=0.5*(xyz2[1]+xyz3[1]);
553     z=0.5*(xyz2[2]+xyz3[2]);
554     }
555     lstpt.push_back(x);
556     lstpt.push_back(y);
557     lstpt.push_back(z);
558     {
559     x=0.5*(xyz2[0]+xyz4[0]);
560     y=0.5*(xyz2[1]+xyz4[1]);
561     z=0.5*(xyz2[2]+xyz4[2]);
562     }
563     lstpt.push_back(x);
564     lstpt.push_back(y);
565     lstpt.push_back(z);
566     {
567     x=0.5*(xyz3[0]+xyz4[0]);
568     y=0.5*(xyz3[1]+xyz4[1]);
569     z=0.5*(xyz3[2]+xyz4[2]);
570     }
571     lstpt.push_back(x);
572     lstpt.push_back(y);
573     lstpt.push_back(z);
574     }