ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_triangle3.cpp
Revision: 638
Committed: Mon Jan 26 21:56:20 2015 UTC (10 years, 7 months ago) by francois
File size: 16492 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_triangle3.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_triangle3.h"
27 #include "fem_maillage.h"
28 #include "fem_noeud.h"
29 #include "mg_element_maillage.h"
30 #include "math.h"
31
32
33 FEM_TRIANGLE3::FEM_TRIANGLE3(unsigned long num,class MG_ELEMENT_MAILLAGE* mai,class FEM_NOEUD** tabnoeud):FEM_ELEMENT2(num,mai),FEM_TEMPLATE_ELEMENT<3>(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 get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
40 }
41
42 FEM_TRIANGLE3::FEM_TRIANGLE3(class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(mai),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
43 {
44 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
45 tab[0]->get_lien_element2()->ajouter(this);
46 tab[1]->get_lien_element2()->ajouter(this);
47 tab[2]->get_lien_element2()->ajouter(this);
48 get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
49 }
50 FEM_TRIANGLE3::FEM_TRIANGLE3(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo,class FEM_NOEUD** tabnoeud):FEM_ELEMENT2(num,topo),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
51 {
52 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
53 tab[0]->get_lien_element2()->ajouter(this);
54 tab[1]->get_lien_element2()->ajouter(this);
55 tab[2]->get_lien_element2()->ajouter(this);
56 get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
57 }
58
59 FEM_TRIANGLE3::FEM_TRIANGLE3(class MG_ELEMENT_TOPOLOGIQUE* topo,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(topo),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
60 {
61 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
62 tab[0]->get_lien_element2()->ajouter(this);
63 tab[1]->get_lien_element2()->ajouter(this);
64 tab[2]->get_lien_element2()->ajouter(this);
65 get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
66 }
67 FEM_TRIANGLE3::FEM_TRIANGLE3(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<3>(tabnoeud)
68 {
69 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
70 tab[0]->get_lien_element2()->ajouter(this);
71 tab[1]->get_lien_element2()->ajouter(this);
72 tab[2]->get_lien_element2()->ajouter(this);
73 get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
74 }
75
76 FEM_TRIANGLE3::FEM_TRIANGLE3(class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(topo,mai),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
77 {
78 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
79 tab[0]->get_lien_element2()->ajouter(this);
80 tab[1]->get_lien_element2()->ajouter(this);
81 tab[2]->get_lien_element2()->ajouter(this);
82 get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
83 }
84 FEM_TRIANGLE3::FEM_TRIANGLE3(FEM_TRIANGLE3& mdd):FEM_ELEMENT2(mdd),FEM_TEMPLATE_ELEMENT<3>(mdd)
85 {
86 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
87 tab[0]->get_lien_element2()->ajouter(this);
88 tab[1]->get_lien_element2()->ajouter(this);
89 tab[2]->get_lien_element2()->ajouter(this);
90 get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
91 }
92 FEM_TRIANGLE3::~FEM_TRIANGLE3()
93 {
94 if (liaison_topologique==NULL) return;
95 if (liaison_topologique->get_dimension()==0) liaison_topologique->get_lien_fem_maillage()->supprimer(this);
96 tab[0]->get_lien_element2()->supprimer(this);
97 tab[1]->get_lien_element2()->supprimer(this);
98 tab[2]->get_lien_element2()->supprimer(this);
99 get_fem_noeudpetitid()->get_lien_petit_element2()->supprimer(this);
100 }
101
102
103 FEM_ELEMENT_MAILLAGE* FEM_TRIANGLE3::dupliquer(FEM_MAILLAGE *femmai,long decalage)
104 {
105 FEM_NOEUD* tabnoeud[3];
106 tabnoeud[0]=femmai->get_fem_noeudid(tab[0]->get_id()+decalage);
107 tabnoeud[1]=femmai->get_fem_noeudid(tab[1]->get_id()+decalage);
108 tabnoeud[2]=femmai->get_fem_noeudid(tab[2]->get_id()+decalage);
109 FEM_TRIANGLE3* tri=new FEM_TRIANGLE3(get_id()+decalage,maillage,tabnoeud);
110 femmai->ajouter_fem_element2(tri);
111 return tri;
112 }
113
114
115
116
117
118 int FEM_TRIANGLE3::get_type_entite(void)
119 {
120 return IDFEM_TRIANGLE3;
121 }
122
123 int FEM_TRIANGLE3::get_dimension(void)
124 {
125 return 2;
126 }
127
128
129
130 int FEM_TRIANGLE3::get_nb_fem_noeud(void)
131 {
132 return FEM_TEMPLATE_ELEMENT<3>::get_nb_fem_noeud();
133 }
134
135 FEM_NOEUD* FEM_TRIANGLE3::get_fem_noeud(int num)
136 {
137 return FEM_TEMPLATE_ELEMENT<3>::get_fem_noeud(num);
138 }
139
140 void FEM_TRIANGLE3::change_noeud(int num,FEM_NOEUD* noeud)
141 {
142 FEM_TEMPLATE_ELEMENT<3>::change_noeud(num,noeud);
143 }
144
145 BOITE_3D& FEM_TRIANGLE3::get_boite_3D(void)
146 {
147 return FEM_TEMPLATE_ELEMENT<3>::get_boite_3D();
148 }
149
150
151 void FEM_TRIANGLE3::enregistrer(std::ostream& o)
152 {
153 if (maillage!=NULL)
154 if (get_lien_topologie()!=NULL) o << "%" << get_id() << "=FEM_TRIANGLE3($"<< get_lien_topologie()->get_id() << ",$" << maillage->get_id() << ",$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<< ");" << std::endl;
155 else o << "%" << get_id() << "=FEM_TRIANGLE3(NULL,$" << maillage->get_id() << ",$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<< ");" << std::endl;
156 else
157 if (get_lien_topologie()!=NULL) o << "%" << get_id() << "=FEM_TRIANGLE3($"<< get_lien_topologie()->get_id() << ",NULL,$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<< ");" << std::endl;
158 else o << "%" << get_id() << "=FEM_TRIANGLE3(NULL,NULL,$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<< ");" << std::endl;
159
160 }
161
162 int FEM_TRIANGLE3::get_nb_pt_gauss(int degre)
163 {
164 if (degre<2) return 1;
165 if (degre<3) return 3;
166 if (degre<4) return 4;
167 if (degre<5) return 6;
168 if (degre<6) return 7;
169 return 0;
170 }
171
172 void FEM_TRIANGLE3::get_pt_gauss(int degre,int num,double &w,double *uv)
173 {
174 if (degre<2)
175 {
176 if (num==0) {w=0.5;uv[0]=0.333333333333333;uv[1]=0.333333333333333;}
177 return;
178 }
179 if (degre<3)
180 {
181 if (num==0) {w=0.166666666666667;uv[0]=0.166666666666667;uv[1]=0.166666666666667;}
182 if (num==1) {w=0.166666666666667;uv[0]=0.666666666666667;uv[1]=0.166666666666667;}
183 if (num==2) {w=0.166666666666667;uv[0]=0.166666666666667;uv[1]=0.666666666666667;}
184 return;
185 }
186 if (degre<4)
187 {
188 if (num==0) {w=-0.28125;uv[0]=0.333333333333333;uv[1]=0.333333333333333;}
189 if (num==1) {w=0.260416666666667;uv[0]=0.2;uv[1]=0.2;}
190 if (num==2) {w=0.260416666666667;uv[0]=0.6;uv[1]=0.2;}
191 if (num==3) {w=0.260416666666667;uv[0]=0.2;uv[1]=0.6;}
192 return;
193 }
194 if (degre<5)
195 {
196 if (num==0) {w=0.111690794839005;uv[0]=0.445948490915965;uv[1]=0.445948490915965;}
197 if (num==1) {w=0.111690794839005;uv[0]=0.108103018168070;uv[1]=0.445948490915965;}
198 if (num==2) {w=0.111690794839005;uv[0]=0.445948490915965;uv[1]=0.108103018168070;}
199 if (num==3) {w=0.054975871827661;uv[0]=0.091576213509771;uv[1]=0.091576213509771;}
200 if (num==4) {w=0.054975871827661;uv[0]=0.816847572980458;uv[1]=0.091576213509771;}
201 if (num==5) {w=0.054975871827661;uv[0]=0.091576213509771;uv[1]=0.816847572980458;}
202 return;
203 }
204 if (degre<6)
205 {
206 if (num==0) {w=0.1125;uv[0]=0.333333333333333;uv[1]=0.333333333333333;}
207 if (num==1) {w=0.066197076394253;uv[0]=0.470142064105115;uv[1]=0.470142064105115;}
208 if (num==2) {w=0.066197076394253;uv[0]=0.059715871789770;uv[1]=0.470142064105115;}
209 if (num==3) {w=0.066197076394253;uv[0]=0.470142064105115;uv[1]=0.059715871789770;}
210 if (num==4) {w=0.062969590272413;uv[0]=0.101286507323456;uv[1]=0.101286507323456;}
211 if (num==5) {w=0.062969590272413;uv[0]=0.797426985353088;uv[1]=0.101286507323456;}
212 if (num==6) {w=0.062969590272413;uv[0]=0.101286507323456;uv[1]=0.797426985353088;}
213 return;
214 }
215 }
216
217 int FEM_TRIANGLE3::get_nb_fonction_interpolation(void)
218 {
219 return 3;
220 }
221 int FEM_TRIANGLE3::get_degremax_fonction_interpolation(void)
222 {
223 return 1;
224 }
225 double FEM_TRIANGLE3::get_fonction_interpolation(int num,double *uv)
226 {
227 double val;
228 switch (num)
229 {
230 case 1:
231 val=1-uv[0]-uv[1];
232 break;
233 case 2:
234 val=uv[0];
235 break;
236 case 3:
237 val=uv[1];
238 break;
239
240 }
241 return val;
242 }
243
244 double FEM_TRIANGLE3::get_fonction_derive_interpolation(int num,int num_variable,double *uv)
245 {
246 double val;
247
248 switch (num)
249 {
250 case 1:
251 switch (num_variable)
252 {
253 case 1:
254 val=-1;
255 break;
256 case 2:
257 val=-1;
258 break;
259 } break;
260 case 2:
261 switch (num_variable)
262 {
263 case 1:
264 val=1;
265 break;
266 case 2:
267 val=0.;
268 break;
269 }break;
270 case 3:
271 switch (num_variable)
272 {
273 case 1:
274 val=0.;
275 break;
276 case 2:
277 val=1.;
278 break;
279 }break;
280
281 }
282 return val;
283 }
284
285 double FEM_TRIANGLE3::get_jacobien(double* jac,double *uv,int& li,int& col,double unite)
286 {
287 // Il faut exprimer les coordonnées x,y,z du repère 3D initial Ri dans un repère 2D final Rf de coordonnées
288 // x',y',z' où z' est nul afin de pouvoir utiliser les fonctions d'interpolation de l'élément de référence 2D
289 FEM_NOEUD* n1=get_fem_noeud(0);
290 FEM_NOEUD* n2=get_fem_noeud(1);
291 FEM_NOEUD* n3=get_fem_noeud(2);
292
293 double xyzn1[3];
294 double xyzn2[3];
295 double xyzn3[3];
296
297 xyzn1[0]=n1->get_x();
298 xyzn1[1]=n1->get_y();
299 xyzn1[2]=n1->get_z();
300
301 xyzn2[0]=n2->get_x();
302 xyzn2[1]=n2->get_y();
303 xyzn2[2]=n2->get_z();
304
305 xyzn3[0]=n3->get_x();
306 xyzn3[1]=n3->get_y();
307 xyzn3[2]=n3->get_z();
308
309 // Position des noeuds dans le repère initial (i)
310 OT_VECTEUR_3D n1_i(xyzn1);
311 OT_VECTEUR_3D n2_i(xyzn2);
312 OT_VECTEUR_3D n3_i(xyzn3);
313 OT_MATRICE_3D pos_i(n1_i,n2_i,n3_i);
314
315 // Le nouveau repère est défini par les vecteurs vec1 et vec2 et son origine est le noeud n1
316 OT_VECTEUR_3D vec1(xyzn1,xyzn2);
317 OT_VECTEUR_3D vec3(xyzn1,xyzn3);
318 vec1.norme();
319 OT_VECTEUR_3D n=vec1&vec3; // Normale au triangle au noeud n1
320 n.norme();
321 OT_VECTEUR_3D vec2=n&vec1; // Vecteur normal à vec1 dans le plan du triangle
322 vec2.norme();
323 OT_MATRICE_3D pif(vec1,vec2,n); // Matrice de passage du repère Ri au repère Rf
324
325 // Position des noeuds dans le repère final (f)
326 OT_MATRICE_3D pif_inverse=pif.inverse();
327 OT_MATRICE_3D pos_f=pif_inverse*pos_i;
328 OT_VECTEUR_3D vecn1=pos_f.get_vecteur1();
329 pos_f.change_vecteur1(pos_f.get_vecteur1()-vecn1); // Étant donné que le noeud 1 est l'origine du repère final, on soustrait ses coordonnées à chaque noeud
330 pos_f.change_vecteur2(pos_f.get_vecteur2()-vecn1);
331 pos_f.change_vecteur3(pos_f.get_vecteur3()-vecn1);
332
333 // Jacobien dans le repère final
334 OT_VECTEUR_3D dNi_de(get_fonction_derive_interpolation(1,1,uv),get_fonction_derive_interpolation(2,1,uv),get_fonction_derive_interpolation(3,1,uv));
335 OT_VECTEUR_3D dNi_dn(get_fonction_derive_interpolation(1,2,uv),get_fonction_derive_interpolation(2,2,uv),get_fonction_derive_interpolation(3,2,uv));
336 OT_VECTEUR_3D dNi_dZ(0.,0.,0.);
337
338 OT_MATRICE_3D dNi(dNi_de,dNi_dn,dNi_dZ);
339 OT_MATRICE_3D dNi_transpose=dNi.transpose();
340 OT_MATRICE_3D pos_f_transpose=pos_f.transpose();
341 OT_MATRICE_3D J_f=dNi_transpose*pos_f_transpose;
342
343 // Jacobien dans le repère initial
344 OT_MATRICE_3D J_i=pif*J_f;
345
346 jac[0]=J_i(0,0);
347 jac[1]=J_i(1,0);
348 jac[2]=J_i(2,0);
349
350 jac[3]=J_i(0,1);
351 jac[4]=J_i(1,1);
352 jac[5]=J_i(2,1);
353
354 jac[6]=0.;
355 jac[7]=0.;
356 jac[8]=0.;
357 }
358
359
360 bool FEM_TRIANGLE3::valide_parametre_EF(double* uvw)
361 {
362 if (uvw[0]>=-1e-10)
363 if (uvw[1]>=-1e-10)
364 if (uvw[0]+uvw[1]<=1.+1e-10)
365 return true;
366 return false;
367 }
368
369 void FEM_TRIANGLE3::get_inverse_jacob(double* j,double *uv,double unite)
370 {
371 // Il faut exprimer les coordonnées x,y,z du repère 3D initial Ri dans un repère 2D final Rf de coordonnées
372 // x',y',z' où z' est nul afin de pouvoir utiliser les fonctions d'interpolation de l'élément de référence
373 FEM_NOEUD* n1=get_fem_noeud(0);
374 FEM_NOEUD* n2=get_fem_noeud(1);
375 FEM_NOEUD* n3=get_fem_noeud(2);
376
377 double xyzn1[3];
378 double xyzn2[3];
379 double xyzn3[3];
380
381 xyzn1[0]=n1->get_x();
382 xyzn1[1]=n1->get_y();
383 xyzn1[2]=n1->get_z();
384
385 xyzn2[0]=n2->get_x();
386 xyzn2[1]=n2->get_y();
387 xyzn2[2]=n2->get_z();
388
389 xyzn3[0]=n3->get_x();
390 xyzn3[1]=n3->get_y();
391 xyzn3[2]=n3->get_z();
392
393 // Position des noeuds dans le repère initial (i)
394 OT_VECTEUR_3D n1_i(xyzn1);
395 OT_VECTEUR_3D n2_i(xyzn2);
396 OT_VECTEUR_3D n3_i(xyzn3);
397 OT_MATRICE_3D pos_i(n1_i,n2_i,n3_i);
398
399 // Le nouveau repère est défini par les vecteurs vec1 et vec2 et son origine est le noeud n1
400 OT_VECTEUR_3D vec1(xyzn1,xyzn2);
401 OT_VECTEUR_3D vec3(xyzn1,xyzn3);
402 vec1.norme();
403 OT_VECTEUR_3D n=vec1&vec3; // Normale au triangle au noeud n1
404 n.norme();
405 OT_VECTEUR_3D vec2=n&vec1; // Vecteur normal à vec1 dans le plan du triangle
406 vec2.norme();
407 OT_MATRICE_3D pif(vec1,vec2,n); // Matrice de passage du repère Ri au repère Rf
408
409 // Position des noeuds dans le repère final (f)
410 OT_MATRICE_3D pif_inverse=pif.inverse();
411 OT_MATRICE_3D pos_f=pif_inverse*pos_i;
412 OT_VECTEUR_3D vecn1=pos_f.get_vecteur1();
413 pos_f.change_vecteur1(pos_f.get_vecteur1()-vecn1); // Étant donné que le noeud 1 est l'origine du repère final, on soustrait ses coordonnées à chaque noeud
414 pos_f.change_vecteur2(pos_f.get_vecteur2()-vecn1);
415 pos_f.change_vecteur3(pos_f.get_vecteur3()-vecn1);
416
417 // Déterminant du jacobien 2D dans le repère final
418 OT_VECTEUR_3D dNi_de(get_fonction_derive_interpolation(1,1,uv),get_fonction_derive_interpolation(2,1,uv),get_fonction_derive_interpolation(3,1,uv));
419 OT_VECTEUR_3D dNi_dn(get_fonction_derive_interpolation(1,2,uv),get_fonction_derive_interpolation(2,2,uv),get_fonction_derive_interpolation(3,2,uv));
420 OT_VECTEUR_3D dNi_dZ(0.,0.,0.);
421
422 OT_MATRICE_3D dNi(dNi_de,dNi_dn,dNi_dZ);
423 OT_MATRICE_3D dNi_transpose=dNi.transpose();
424 OT_MATRICE_3D pos_f_transpose=pos_f.transpose();
425 OT_MATRICE_3D J=dNi_transpose*pos_f_transpose;
426
427 double J11=J(0,0);
428 double J21=J(1,0);
429 double J12=J(0,1);
430 double J22=J(1,1);
431 double det_J=J11*J22-J12*J21;
432
433 // Jacobien inverse 3D dans le repère final
434 OT_VECTEUR_3D col1_jf_3d(J22/det_J,-J21/det_J,0.);
435 OT_VECTEUR_3D col2_jf_3d(-J12/det_J,J11/det_J,0.);
436 OT_VECTEUR_3D col3_jf_3d(0.,0.,0.);
437 OT_MATRICE_3D j_f_3d(col1_jf_3d,col2_jf_3d,col3_jf_3d);
438
439 // Jacobien inverse 3D dans le repère initial
440 OT_MATRICE_3D j_i=pif*j_f_3d;
441
442 j[0]=j_i(0,0);
443 j[1]=j_i(1,0);
444 j[2]=j_i(2,0);
445
446 j[3]=j_i(0,1);
447 j[4]=j_i(1,1);
448 j[5]=j_i(2,1);
449
450 j[6]=j_i(0,2);
451 j[7]=j_i(1,2);
452 j[8]=j_i(2,2);
453 }