ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_triangle6.cpp
Revision: 638
Committed: Mon Jan 26 21:56:20 2015 UTC (10 years, 7 months ago) by francois
File size: 16683 byte(s)
Log Message:
ajout d'une méthode qui renvoit le degre max des fonctions d'interpolation et changement de nom de la fonction qui renvoie le nombre de fonction d'interpolation (ajout d'un get pour uniformiser)

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::get_nb_fonction_interpolation(void)
258 {
259 return 6;
260 }
261 int FEM_TRIANGLE6::get_degremax_fonction_interpolation(void)
262 {
263 return 2;
264 }
265 double FEM_TRIANGLE6::get_fonction_interpolation(int num,double *uv)
266 {
267 double val;
268 switch (num)
269 {
270 case 1:
271 val=(-1+uv[0]+uv[1])*(-1+2*uv[0]+2*uv[1]);
272 break;
273 case 2:
274 val=4.*uv[0]*(1-uv[0]-uv[1]);
275 break;
276 case 3:
277 val=-uv[0]*(1.-2.*uv[0]);
278 break;
279 case 4:
280 val=4.*uv[0]*uv[1];
281 break;
282 case 5:
283 val=-uv[1]*(1-2*uv[1]);
284 break;
285 case 6:
286 val=4.*uv[1]*(1-uv[0]-uv[1]);
287 break;
288
289 }
290 return val;
291 }
292
293 double FEM_TRIANGLE6::get_fonction_derive_interpolation(int num,int num_variable,double *uv)
294 {
295 double val;
296
297 switch (num)
298 {
299 case 1:
300 switch (num_variable)
301 {
302 case 1:
303 val=-3.+4*uv[0]+4*uv[1];
304 break;
305 case 2:
306 val=-3.+4*uv[0]+4*uv[1];
307 break;
308 } break;
309 case 2:
310 switch (num_variable)
311 {
312 case 1:
313 val=4.-8*uv[0]-4*uv[1];
314 break;
315 case 2:
316 val=-4.*uv[0];
317 break;
318 }break;
319 case 3:
320 switch (num_variable)
321 {
322 case 1:
323 val=-1.+4*uv[0];
324 break;
325 case 2:
326 val=0.;
327 break;
328 }break;
329 case 4:
330 switch (num_variable)
331 {
332 case 1:
333 val=4*uv[1];
334 break;
335 case 2:
336 val=4*uv[0];
337 break;
338 }break;
339 case 5:
340 switch (num_variable)
341 {
342 case 1:
343 val=0.;
344 break;
345 case 2:
346 val=-1+4.*uv[1];
347 break;
348 }break;
349 case 6:
350 switch (num_variable)
351 {
352 case 1:
353 val=-4*uv[1];
354 break;
355 case 2:
356 val=4-4*uv[0]-8*uv[1];
357 break;
358 }break;
359
360 }
361 return val;
362 }
363
364
365 double FEM_TRIANGLE6::get_jacobien(double* jac,double *uv,int& li,int& col,double unite)
366 {
367 FEM_NOEUD* n1=get_fem_noeud(0);
368 FEM_NOEUD* n2=get_fem_noeud(1);
369 FEM_NOEUD* n3=get_fem_noeud(2);
370 FEM_NOEUD* n4=get_fem_noeud(3);
371 FEM_NOEUD* n5=get_fem_noeud(4);
372 FEM_NOEUD* n6=get_fem_noeud(5);
373
374 double xyzn1[3],xyzn2[3],xyzn3[3],xyzn4[3],xyzn5[3],xyzn6[3];
375
376 xyzn1[0]=n1->get_x();
377 xyzn1[1]=n1->get_y();
378 xyzn1[2]=n1->get_z();
379
380 xyzn2[0]=n2->get_x();
381 xyzn2[1]=n2->get_y();
382 xyzn2[2]=n2->get_z();
383
384 xyzn3[0]=n3->get_x();
385 xyzn3[1]=n3->get_y();
386 xyzn3[2]=n3->get_z();
387
388 xyzn4[0]=n4->get_x();
389 xyzn4[1]=n4->get_y();
390 xyzn4[2]=n4->get_z();
391
392 xyzn5[0]=n5->get_x();
393 xyzn5[1]=n5->get_y();
394 xyzn5[2]=n5->get_z();
395
396 xyzn6[0]=n6->get_x();
397 xyzn6[1]=n6->get_y();
398 xyzn6[2]=n6->get_z();
399
400 OT_VECTEUR_3D vec1(xyzn1,xyzn3);
401 OT_VECTEUR_3D vec3(xyzn1,xyzn5);
402 OT_VECTEUR_3D n=vec1&vec3;
403 n.norme();
404
405 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];
406 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];
407 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];
408
409 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];
410 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];
411 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];
412
413 jac[6]=n[0];
414 jac[7]=n[1];
415 jac[8]=n[2];
416
417 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]);
418
419 return DET;
420
421 }
422
423 void FEM_TRIANGLE6::get_inverse_jacob(double* j,double *uv,double unite)
424 {
425 double jac[9];
426 int lig,col;
427 double detj=this->get_jacobien(jac,uv,lig,col,unite);
428
429 j[0] = (jac[4]*jac[8]-jac[5]*jac[7])/(detj);
430 j[1] = -(jac[1]*jac[8]-jac[2]*jac[7])/(detj);
431 j[2] =-(-jac[1]*jac[5]+jac[2]*jac[4])/(detj);
432
433 j[3] = -(jac[1]*jac[8]-jac[5]*jac[6])/(detj);
434 j[4] = (jac[0]*jac[8]-jac[2]*jac[6])/(detj);
435 j[5] = -(jac[0]*jac[5]-jac[2]*jac[3])/(detj);
436
437 j[6] =-(-jac[3]*jac[7]+jac[4]*jac[6])/(detj);
438 j[7] = -(jac[0]*jac[7]-jac[1]*jac[6])/(detj);
439 j[8] = (jac[0]*jac[4]-jac[1]*jac[4])/(detj);
440
441 }
442