ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/dly_tetra.cpp
Revision: 469
Committed: Wed Nov 27 20:02:22 2013 UTC (11 years, 5 months ago) by francois
Original Path: magic/lib/mailleur_auto/src/dly_tetra.cpp
File size: 15633 byte(s)
Log Message:
ajout d'un generateur de carte de taille constante pour les nouvelles cartes et decoupage de la librairie mailleur en 2 pour cause d'appel cyclique

File Contents

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