ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_triangle6.cpp
Revision: 637
Committed: Mon Jan 26 21:26:37 2015 UTC (10 years, 3 months ago) by francois
File size: 16598 byte(s)
Log Message:
entree des points de gauss pour les différents éléments finis depuis reference JCC

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     // fem_triangle6.cpp
16     //
17     //------------------------------------------------------------
18     //------------------------------------------------------------
19     // COPYRIGHT 2000
20     // Version du 02/03/2006 � 11H22
21     //------------------------------------------------------------
22     //------------------------------------------------------------
23    
24    
25     #include "gestionversion.h"
26     #include "fem_triangle6.h"
27     #include "fem_noeud.h"
28     #include "mg_element_maillage.h"
29     #include "fem_maillage.h"
30 francois 612 #include "math.h"
31     #include "ot_quadrature_gauss.h"
32 francois 283
33 francois 309 FEM_TRIANGLE6::FEM_TRIANGLE6(unsigned long num,class MG_ELEMENT_MAILLAGE* mai,class FEM_NOEUD** tabnoeud):FEM_ELEMENT2(num,mai),FEM_TEMPLATE_ELEMENT<6>(tabnoeud)
34 francois 283 {
35     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
36 francois 309 tab[0]->get_lien_element2()->ajouter(this);
37     tab[1]->get_lien_element2()->ajouter(this);
38     tab[2]->get_lien_element2()->ajouter(this);
39     tab[3]->get_lien_element2()->ajouter(this);
40     tab[4]->get_lien_element2()->ajouter(this);
41     tab[5]->get_lien_element2()->ajouter(this);
42     get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
43 francois 283 }
44    
45 francois 309 FEM_TRIANGLE6::FEM_TRIANGLE6(class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(mai),FEM_TEMPLATE_ELEMENT<6>(tabnoeud)
46 francois 283 {
47     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
48 francois 309 tab[0]->get_lien_element2()->ajouter(this);
49     tab[1]->get_lien_element2()->ajouter(this);
50     tab[2]->get_lien_element2()->ajouter(this);
51     tab[3]->get_lien_element2()->ajouter(this);
52     tab[4]->get_lien_element2()->ajouter(this);
53     tab[5]->get_lien_element2()->ajouter(this);
54     get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
55 francois 283 }
56    
57 francois 378 FEM_TRIANGLE6::FEM_TRIANGLE6(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo,class FEM_NOEUD** tabnoeud):FEM_ELEMENT2(num,topo),FEM_TEMPLATE_ELEMENT<6>(tabnoeud)
58     {
59     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
60     tab[0]->get_lien_element2()->ajouter(this);
61     tab[1]->get_lien_element2()->ajouter(this);
62     tab[2]->get_lien_element2()->ajouter(this);
63     tab[3]->get_lien_element2()->ajouter(this);
64     tab[4]->get_lien_element2()->ajouter(this);
65     tab[5]->get_lien_element2()->ajouter(this);
66     get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
67     }
68    
69     FEM_TRIANGLE6::FEM_TRIANGLE6(class MG_ELEMENT_TOPOLOGIQUE* topo,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(topo),FEM_TEMPLATE_ELEMENT<6>(tabnoeud)
70     {
71     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
72     tab[0]->get_lien_element2()->ajouter(this);
73     tab[1]->get_lien_element2()->ajouter(this);
74     tab[2]->get_lien_element2()->ajouter(this);
75     tab[3]->get_lien_element2()->ajouter(this);
76     tab[4]->get_lien_element2()->ajouter(this);
77     tab[5]->get_lien_element2()->ajouter(this);
78     get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
79     }
80     FEM_TRIANGLE6::FEM_TRIANGLE6(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai,class FEM_NOEUD** tabnoeud):FEM_ELEMENT2(num,topo,mai),FEM_TEMPLATE_ELEMENT<6>(tabnoeud)
81     {
82     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
83     tab[0]->get_lien_element2()->ajouter(this);
84     tab[1]->get_lien_element2()->ajouter(this);
85     tab[2]->get_lien_element2()->ajouter(this);
86     tab[3]->get_lien_element2()->ajouter(this);
87     tab[4]->get_lien_element2()->ajouter(this);
88     tab[5]->get_lien_element2()->ajouter(this);
89     get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
90     }
91    
92     FEM_TRIANGLE6::FEM_TRIANGLE6(class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(topo,mai),FEM_TEMPLATE_ELEMENT<6>(tabnoeud)
93     {
94     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
95     tab[0]->get_lien_element2()->ajouter(this);
96     tab[1]->get_lien_element2()->ajouter(this);
97     tab[2]->get_lien_element2()->ajouter(this);
98     tab[3]->get_lien_element2()->ajouter(this);
99     tab[4]->get_lien_element2()->ajouter(this);
100     tab[5]->get_lien_element2()->ajouter(this);
101     get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
102     }
103    
104    
105    
106    
107    
108    
109    
110    
111 francois 309 FEM_TRIANGLE6::FEM_TRIANGLE6(FEM_TRIANGLE6& mdd):FEM_ELEMENT2(mdd),FEM_TEMPLATE_ELEMENT<6>(mdd)
112 francois 283 {
113     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
114 francois 309 tab[0]->get_lien_element2()->ajouter(this);
115     tab[1]->get_lien_element2()->ajouter(this);
116     tab[2]->get_lien_element2()->ajouter(this);
117     tab[3]->get_lien_element2()->ajouter(this);
118     tab[4]->get_lien_element2()->ajouter(this);
119     tab[5]->get_lien_element2()->ajouter(this);
120     get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
121 francois 283 }
122     FEM_TRIANGLE6::~FEM_TRIANGLE6()
123     {
124     if (liaison_topologique==NULL) return;
125     if (liaison_topologique->get_dimension()==0) liaison_topologique->get_lien_fem_maillage()->supprimer(this);
126 francois 309 tab[0]->get_lien_element2()->supprimer(this);
127     tab[1]->get_lien_element2()->supprimer(this);
128     tab[2]->get_lien_element2()->supprimer(this);
129     tab[3]->get_lien_element2()->supprimer(this);
130     tab[4]->get_lien_element2()->supprimer(this);
131     tab[5]->get_lien_element2()->supprimer(this);
132     get_fem_noeudpetitid()->get_lien_petit_element2()->supprimer(this);
133 francois 283 }
134    
135     FEM_ELEMENT_MAILLAGE* FEM_TRIANGLE6::dupliquer(FEM_MAILLAGE *femmai,long decalage)
136     {
137     FEM_NOEUD* tabnoeud[6];
138     tabnoeud[0]=femmai->get_fem_noeudid(tab[0]->get_id()+decalage);
139     tabnoeud[1]=femmai->get_fem_noeudid(tab[1]->get_id()+decalage);
140     tabnoeud[2]=femmai->get_fem_noeudid(tab[2]->get_id()+decalage);
141     tabnoeud[3]=femmai->get_fem_noeudid(tab[3]->get_id()+decalage);
142     tabnoeud[4]=femmai->get_fem_noeudid(tab[4]->get_id()+decalage);
143     tabnoeud[5]=femmai->get_fem_noeudid(tab[5]->get_id()+decalage);
144     FEM_TRIANGLE6* tri=new FEM_TRIANGLE6(get_id()+decalage,maillage,tabnoeud);
145 francois 309 femmai->ajouter_fem_element2(tri);
146 francois 283 return tri;
147     }
148    
149 francois 635 bool FEM_TRIANGLE6::valide_parametre_EF(double* uvw)
150     {
151     if (uvw[0]>=-1e-10)
152     if (uvw[1]>=-1e-10)
153     if (uvw[0]+uvw[1]<=1.+1e-10)
154     return true;
155     return false;
156     }
157 francois 283
158    
159     int FEM_TRIANGLE6::get_type_entite(void)
160     {
161     return IDFEM_TRIANGLE6;
162     }
163    
164     int FEM_TRIANGLE6::get_dimension(void)
165     {
166     return 2;
167     }
168    
169    
170     int FEM_TRIANGLE6::get_nb_fem_noeud(void)
171     {
172     return FEM_TEMPLATE_ELEMENT<6>::get_nb_fem_noeud();
173     }
174    
175     FEM_NOEUD* FEM_TRIANGLE6::get_fem_noeud(int num)
176     {
177     return FEM_TEMPLATE_ELEMENT<6>::get_fem_noeud(num);
178     }
179    
180     void FEM_TRIANGLE6::change_noeud(int num,FEM_NOEUD* noeud)
181     {
182     FEM_TEMPLATE_ELEMENT<6>::change_noeud(num,noeud);
183     }
184    
185     BOITE_3D& FEM_TRIANGLE6::get_boite_3D(void)
186     {
187     return FEM_TEMPLATE_ELEMENT<6>::get_boite_3D();
188     }
189    
190    
191     void FEM_TRIANGLE6::enregistrer(std::ostream& o)
192     {
193 francois 378 if (maillage!=NULL)
194     if (get_lien_topologie()!=NULL) o << "%" << get_id() << "=FEM_TRIANGLE6($"<< get_lien_topologie()->get_id() << ",$" << maillage->get_id() << ",$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id() << ",$" << tab[3]->get_id() << ",$" << tab[4]->get_id() << ",$" << tab[5]->get_id() << ");" << std::endl;
195     else o << "%" << get_id() << "=FEM_TRIANGLE6(NULL,$" << maillage->get_id() << ",$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id() << ",$" << tab[3]->get_id() << ",$" << tab[4]->get_id() << ",$" << tab[5]->get_id() << ");" << std::endl;
196     else
197     if (get_lien_topologie()!=NULL) o << "%" << get_id() << "=FEM_TRIANGLE6($"<< get_lien_topologie()->get_id() << ",NULL,$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id() << ",$" << tab[3]->get_id() << ",$" << tab[4]->get_id() << ",$" << tab[5]->get_id() << ");" << std::endl;
198     else o << "%" << get_id() << "=FEM_TRIANGLE6(NULL,NULL,$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id() << ",$" << tab[3]->get_id() << ",$" << tab[4]->get_id() << ",$" << tab[5]->get_id() << ");" << std::endl;
199    
200 francois 283 }
201 francois 310
202 francois 637 int FEM_TRIANGLE6::get_nb_pt_gauss(int degre)
203     {
204     if (degre<2) return 1;
205     if (degre<3) return 3;
206     if (degre<4) return 4;
207     if (degre<5) return 6;
208     if (degre<6) return 7;
209     return 0;
210     }
211 francois 310
212 francois 637 void FEM_TRIANGLE6::get_pt_gauss(int degre,int num,double &w,double *uv)
213     {
214     if (degre<2)
215     {
216     if (num==0) {w=0.5;uv[0]=0.333333333333333;uv[1]=0.333333333333333;}
217     return;
218     }
219     if (degre<3)
220     {
221     if (num==0) {w=0.166666666666667;uv[0]=0.166666666666667;uv[1]=0.166666666666667;}
222     if (num==1) {w=0.166666666666667;uv[0]=0.666666666666667;uv[1]=0.166666666666667;}
223     if (num==2) {w=0.166666666666667;uv[0]=0.166666666666667;uv[1]=0.666666666666667;}
224     return;
225     }
226     if (degre<4)
227     {
228     if (num==0) {w=-0.28125;uv[0]=0.333333333333333;uv[1]=0.333333333333333;}
229     if (num==1) {w=0.260416666666667;uv[0]=0.2;uv[1]=0.2;}
230     if (num==2) {w=0.260416666666667;uv[0]=0.6;uv[1]=0.2;}
231     if (num==3) {w=0.260416666666667;uv[0]=0.2;uv[1]=0.6;}
232     return;
233     }
234     if (degre<5)
235     {
236     if (num==0) {w=0.111690794839005;uv[0]=0.445948490915965;uv[1]=0.445948490915965;}
237     if (num==1) {w=0.111690794839005;uv[0]=0.108103018168070;uv[1]=0.445948490915965;}
238     if (num==2) {w=0.111690794839005;uv[0]=0.445948490915965;uv[1]=0.108103018168070;}
239     if (num==3) {w=0.054975871827661;uv[0]=0.091576213509771;uv[1]=0.091576213509771;}
240     if (num==4) {w=0.054975871827661;uv[0]=0.816847572980458;uv[1]=0.091576213509771;}
241     if (num==5) {w=0.054975871827661;uv[0]=0.091576213509771;uv[1]=0.816847572980458;}
242     return;
243     }
244     if (degre<6)
245     {
246     if (num==0) {w=0.1125;uv[0]=0.333333333333333;uv[1]=0.333333333333333;}
247     if (num==1) {w=0.066197076394253;uv[0]=0.470142064105115;uv[1]=0.470142064105115;}
248     if (num==2) {w=0.066197076394253;uv[0]=0.059715871789770;uv[1]=0.470142064105115;}
249     if (num==3) {w=0.066197076394253;uv[0]=0.470142064105115;uv[1]=0.059715871789770;}
250     if (num==4) {w=0.062969590272413;uv[0]=0.101286507323456;uv[1]=0.101286507323456;}
251     if (num==5) {w=0.062969590272413;uv[0]=0.797426985353088;uv[1]=0.101286507323456;}
252     if (num==6) {w=0.062969590272413;uv[0]=0.101286507323456;uv[1]=0.797426985353088;}
253     return;
254     }
255     }
256    
257 francois 310 int FEM_TRIANGLE6::nb_fonction_interpolation(void)
258     {
259     return 6;
260     }
261    
262     double FEM_TRIANGLE6::get_fonction_interpolation(int num,double *uv)
263     {
264     double val;
265     switch (num)
266     {
267     case 1:
268     val=(-1+uv[0]+uv[1])*(-1+2*uv[0]+2*uv[1]);
269     break;
270     case 2:
271     val=4.*uv[0]*(1-uv[0]-uv[1]);
272     break;
273     case 3:
274     val=-uv[0]*(1.-2.*uv[0]);
275     break;
276     case 4:
277     val=4.*uv[0]*uv[1];
278     break;
279     case 5:
280     val=-uv[1]*(1-2*uv[1]);
281     break;
282     case 6:
283     val=4.*uv[1]*(1-uv[0]-uv[1]);
284     break;
285    
286     }
287     return val;
288     }
289    
290     double FEM_TRIANGLE6::get_fonction_derive_interpolation(int num,int num_variable,double *uv)
291     {
292     double val;
293    
294     switch (num)
295     {
296     case 1:
297     switch (num_variable)
298     {
299     case 1:
300 francois 612 val=-3.+4*uv[0]+4*uv[1];
301 francois 310 break;
302     case 2:
303     val=-3.+4*uv[0]+4*uv[1];
304     break;
305     } break;
306     case 2:
307     switch (num_variable)
308     {
309     case 1:
310     val=4.-8*uv[0]-4*uv[1];
311     break;
312     case 2:
313     val=-4.*uv[0];
314     break;
315     }break;
316     case 3:
317     switch (num_variable)
318     {
319     case 1:
320 francois 612 val=-1.+4*uv[0];
321 francois 310 break;
322     case 2:
323     val=0.;
324     break;
325     }break;
326     case 4:
327     switch (num_variable)
328     {
329     case 1:
330     val=4*uv[1];
331     break;
332     case 2:
333     val=4*uv[0];
334     break;
335     }break;
336     case 5:
337     switch (num_variable)
338     {
339     case 1:
340     val=0.;
341     break;
342     case 2:
343     val=-1+4.*uv[1];
344     break;
345     }break;
346     case 6:
347     switch (num_variable)
348     {
349     case 1:
350     val=-4*uv[1];
351     break;
352     case 2:
353     val=4-4*uv[0]-8*uv[1];
354     break;
355     }break;
356    
357     }
358     return val;
359     }
360    
361 francois 612
362 francois 310 double FEM_TRIANGLE6::get_jacobien(double* jac,double *uv,int& li,int& col,double unite)
363     {
364 francois 612 FEM_NOEUD* n1=get_fem_noeud(0);
365     FEM_NOEUD* n2=get_fem_noeud(1);
366     FEM_NOEUD* n3=get_fem_noeud(2);
367     FEM_NOEUD* n4=get_fem_noeud(3);
368     FEM_NOEUD* n5=get_fem_noeud(4);
369     FEM_NOEUD* n6=get_fem_noeud(5);
370    
371     double xyzn1[3],xyzn2[3],xyzn3[3],xyzn4[3],xyzn5[3],xyzn6[3];
372    
373     xyzn1[0]=n1->get_x();
374     xyzn1[1]=n1->get_y();
375     xyzn1[2]=n1->get_z();
376    
377     xyzn2[0]=n2->get_x();
378     xyzn2[1]=n2->get_y();
379     xyzn2[2]=n2->get_z();
380    
381     xyzn3[0]=n3->get_x();
382     xyzn3[1]=n3->get_y();
383     xyzn3[2]=n3->get_z();
384    
385     xyzn4[0]=n4->get_x();
386     xyzn4[1]=n4->get_y();
387     xyzn4[2]=n4->get_z();
388    
389     xyzn5[0]=n5->get_x();
390     xyzn5[1]=n5->get_y();
391     xyzn5[2]=n5->get_z();
392    
393     xyzn6[0]=n6->get_x();
394     xyzn6[1]=n6->get_y();
395     xyzn6[2]=n6->get_z();
396    
397     OT_VECTEUR_3D vec1(xyzn1,xyzn3);
398     OT_VECTEUR_3D vec3(xyzn1,xyzn5);
399     OT_VECTEUR_3D n=vec1&vec3;
400     n.norme();
401 francois 310
402 francois 612 jac[0]=get_fonction_derive_interpolation(1,1,uv)*xyzn1[0]+get_fonction_derive_interpolation(2,1,uv)*xyzn2[0]+get_fonction_derive_interpolation(3,1,uv)*xyzn3[0]+get_fonction_derive_interpolation(4,1,uv)*xyzn4[0]+get_fonction_derive_interpolation(6,1,uv)*xyzn6[0];
403     jac[1]=get_fonction_derive_interpolation(1,1,uv)*xyzn1[1]+get_fonction_derive_interpolation(2,1,uv)*xyzn2[1]+get_fonction_derive_interpolation(3,1,uv)*xyzn3[1]+get_fonction_derive_interpolation(4,1,uv)*xyzn4[1]+get_fonction_derive_interpolation(6,1,uv)*xyzn6[1];
404     jac[2]=get_fonction_derive_interpolation(1,1,uv)*xyzn1[2]+get_fonction_derive_interpolation(2,1,uv)*xyzn2[2]+get_fonction_derive_interpolation(3,1,uv)*xyzn3[2]+get_fonction_derive_interpolation(4,1,uv)*xyzn4[2]+get_fonction_derive_interpolation(6,1,uv)*xyzn6[2];
405    
406     jac[3]=get_fonction_derive_interpolation(1,2,uv)*xyzn1[0]+get_fonction_derive_interpolation(2,2,uv)*xyzn2[0]+get_fonction_derive_interpolation(4,2,uv)*xyzn4[0]+get_fonction_derive_interpolation(5,2,uv)*xyzn5[0]+get_fonction_derive_interpolation(6,2,uv)*xyzn6[0];
407     jac[4]=get_fonction_derive_interpolation(1,2,uv)*xyzn1[1]+get_fonction_derive_interpolation(2,2,uv)*xyzn2[1]+get_fonction_derive_interpolation(4,2,uv)*xyzn4[1]+get_fonction_derive_interpolation(5,2,uv)*xyzn5[1]+get_fonction_derive_interpolation(6,2,uv)*xyzn6[1];
408     jac[5]=get_fonction_derive_interpolation(1,2,uv)*xyzn1[2]+get_fonction_derive_interpolation(2,2,uv)*xyzn2[2]+get_fonction_derive_interpolation(4,2,uv)*xyzn4[2]+get_fonction_derive_interpolation(5,2,uv)*xyzn5[2]+get_fonction_derive_interpolation(6,2,uv)*xyzn6[2];
409    
410     jac[6]=n[0];
411     jac[7]=n[1];
412     jac[8]=n[2];
413    
414     double DET= jac[0]*(jac[4]*jac[8]-jac[7]*jac[5])+jac[1]*(jac[5]*jac[6]-jac[8]*jac[3])+jac[2]*(jac[3]*jac[7]-jac[6]*jac[4]);
415    
416     return DET;
417    
418 francois 310 }
419    
420     void FEM_TRIANGLE6::get_inverse_jacob(double* j,double *uv,double unite)
421     {
422 francois 612 double jac[9];
423     int lig,col;
424     double detj=this->get_jacobien(jac,uv,lig,col,unite);
425    
426     j[0] = (jac[4]*jac[8]-jac[5]*jac[7])/(detj);
427     j[1] = -(jac[1]*jac[8]-jac[2]*jac[7])/(detj);
428     j[2] =-(-jac[1]*jac[5]+jac[2]*jac[4])/(detj);
429    
430     j[3] = -(jac[1]*jac[8]-jac[5]*jac[6])/(detj);
431     j[4] = (jac[0]*jac[8]-jac[2]*jac[6])/(detj);
432     j[5] = -(jac[0]*jac[5]-jac[2]*jac[3])/(detj);
433    
434     j[6] =-(-jac[3]*jac[7]+jac[4]*jac[6])/(detj);
435     j[7] = -(jac[0]*jac[7]-jac[1]*jac[6])/(detj);
436     j[8] = (jac[0]*jac[4]-jac[1]*jac[4])/(detj);
437    
438     }
439