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, 7 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

# 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 // 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 #include "math.h"
31 #include "ot_quadrature_gauss.h"
32
33 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 {
35 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
36 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 }
44
45 FEM_TRIANGLE6::FEM_TRIANGLE6(class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(mai),FEM_TEMPLATE_ELEMENT<6>(tabnoeud)
46 {
47 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
48 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 }
56
57 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 FEM_TRIANGLE6::FEM_TRIANGLE6(FEM_TRIANGLE6& mdd):FEM_ELEMENT2(mdd),FEM_TEMPLATE_ELEMENT<6>(mdd)
112 {
113 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
114 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 }
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 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 }
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 femmai->ajouter_fem_element2(tri);
146 return tri;
147 }
148
149 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
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 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 }
201
202 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
212 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 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 val=-3.+4*uv[0]+4*uv[1];
301 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 val=-1.+4*uv[0];
321 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
362 double FEM_TRIANGLE6::get_jacobien(double* jac,double *uv,int& li,int& col,double unite)
363 {
364 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
402 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 }
419
420 void FEM_TRIANGLE6::get_inverse_jacob(double* j,double *uv,double unite)
421 {
422 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