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::nb_fonction_interpolation(void) |
163 |
{ |
164 |
return 3; |
165 |
} |
166 |
|
167 |
double FEM_TRIANGLE3::get_fonction_interpolation(int num,double *uv) |
168 |
{ |
169 |
double val; |
170 |
switch (num) |
171 |
{ |
172 |
case 1: |
173 |
val=1-uv[0]-uv[1]; |
174 |
break; |
175 |
case 2: |
176 |
val=uv[0]; |
177 |
break; |
178 |
case 3: |
179 |
val=uv[1]; |
180 |
break; |
181 |
|
182 |
} |
183 |
return val; |
184 |
} |
185 |
|
186 |
double FEM_TRIANGLE3::get_fonction_derive_interpolation(int num,int num_variable,double *uv) |
187 |
{ |
188 |
double val; |
189 |
|
190 |
switch (num) |
191 |
{ |
192 |
case 1: |
193 |
switch (num_variable) |
194 |
{ |
195 |
case 1: |
196 |
val=-1; |
197 |
break; |
198 |
case 2: |
199 |
val=-1; |
200 |
break; |
201 |
} break; |
202 |
case 2: |
203 |
switch (num_variable) |
204 |
{ |
205 |
case 1: |
206 |
val=1; |
207 |
break; |
208 |
case 2: |
209 |
val=0.; |
210 |
break; |
211 |
}break; |
212 |
case 3: |
213 |
switch (num_variable) |
214 |
{ |
215 |
case 1: |
216 |
val=0.; |
217 |
break; |
218 |
case 2: |
219 |
val=1.; |
220 |
break; |
221 |
}break; |
222 |
|
223 |
} |
224 |
return val; |
225 |
} |
226 |
|
227 |
double FEM_TRIANGLE3::get_jacobien(double* jac,double *uv,int& li,int& col,double unite) |
228 |
{ |
229 |
// 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 |
230 |
// x',y',z' où z' est nul afin de pouvoir utiliser les fonctions d'interpolation de l'élément de référence 2D |
231 |
FEM_NOEUD* n1=get_fem_noeud(0); |
232 |
FEM_NOEUD* n2=get_fem_noeud(1); |
233 |
FEM_NOEUD* n3=get_fem_noeud(2); |
234 |
|
235 |
double xyzn1[3]; |
236 |
double xyzn2[3]; |
237 |
double xyzn3[3]; |
238 |
|
239 |
xyzn1[0]=n1->get_x(); |
240 |
xyzn1[1]=n1->get_y(); |
241 |
xyzn1[2]=n1->get_z(); |
242 |
|
243 |
xyzn2[0]=n2->get_x(); |
244 |
xyzn2[1]=n2->get_y(); |
245 |
xyzn2[2]=n2->get_z(); |
246 |
|
247 |
xyzn3[0]=n3->get_x(); |
248 |
xyzn3[1]=n3->get_y(); |
249 |
xyzn3[2]=n3->get_z(); |
250 |
|
251 |
// Position des noeuds dans le repère initial (i) |
252 |
OT_VECTEUR_3D n1_i(xyzn1); |
253 |
OT_VECTEUR_3D n2_i(xyzn2); |
254 |
OT_VECTEUR_3D n3_i(xyzn3); |
255 |
OT_MATRICE_3D pos_i(n1_i,n2_i,n3_i); |
256 |
|
257 |
// Le nouveau repère est défini par les vecteurs vec1 et vec2 et son origine est le noeud n1 |
258 |
OT_VECTEUR_3D vec1(xyzn1,xyzn2); |
259 |
OT_VECTEUR_3D vec3(xyzn1,xyzn3); |
260 |
vec1.norme(); |
261 |
OT_VECTEUR_3D n=vec1&vec3; // Normale au triangle au noeud n1 |
262 |
n.norme(); |
263 |
OT_VECTEUR_3D vec2=n&vec1; // Vecteur normal à vec1 dans le plan du triangle |
264 |
vec2.norme(); |
265 |
OT_MATRICE_3D pif(vec1,vec2,n); // Matrice de passage du repère Ri au repère Rf |
266 |
|
267 |
// Position des noeuds dans le repère final (f) |
268 |
OT_MATRICE_3D pif_inverse=pif.inverse(); |
269 |
OT_MATRICE_3D pos_f=pif_inverse*pos_i; |
270 |
OT_VECTEUR_3D vecn1=pos_f.get_vecteur1(); |
271 |
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 |
272 |
pos_f.change_vecteur2(pos_f.get_vecteur2()-vecn1); |
273 |
pos_f.change_vecteur3(pos_f.get_vecteur3()-vecn1); |
274 |
|
275 |
// Jacobien dans le repère final |
276 |
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)); |
277 |
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)); |
278 |
OT_VECTEUR_3D dNi_dZ(0.,0.,0.); |
279 |
|
280 |
OT_MATRICE_3D dNi(dNi_de,dNi_dn,dNi_dZ); |
281 |
OT_MATRICE_3D dNi_transpose=dNi.transpose(); |
282 |
OT_MATRICE_3D pos_f_transpose=pos_f.transpose(); |
283 |
OT_MATRICE_3D J_f=dNi_transpose*pos_f_transpose; |
284 |
|
285 |
// Jacobien dans le repère initial |
286 |
OT_MATRICE_3D J_i=pif*J_f; |
287 |
|
288 |
jac[0]=J_i(0,0); |
289 |
jac[1]=J_i(1,0); |
290 |
jac[2]=J_i(2,0); |
291 |
|
292 |
jac[3]=J_i(0,1); |
293 |
jac[4]=J_i(1,1); |
294 |
jac[5]=J_i(2,1); |
295 |
|
296 |
jac[6]=0.; |
297 |
jac[7]=0.; |
298 |
jac[8]=0.; |
299 |
} |
300 |
|
301 |
void FEM_TRIANGLE3::get_inverse_jacob(double* j,double *uv,double unite) |
302 |
{ |
303 |
// 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 |
304 |
// x',y',z' où z' est nul afin de pouvoir utiliser les fonctions d'interpolation de l'élément de référence |
305 |
FEM_NOEUD* n1=get_fem_noeud(0); |
306 |
FEM_NOEUD* n2=get_fem_noeud(1); |
307 |
FEM_NOEUD* n3=get_fem_noeud(2); |
308 |
|
309 |
double xyzn1[3]; |
310 |
double xyzn2[3]; |
311 |
double xyzn3[3]; |
312 |
|
313 |
xyzn1[0]=n1->get_x(); |
314 |
xyzn1[1]=n1->get_y(); |
315 |
xyzn1[2]=n1->get_z(); |
316 |
|
317 |
xyzn2[0]=n2->get_x(); |
318 |
xyzn2[1]=n2->get_y(); |
319 |
xyzn2[2]=n2->get_z(); |
320 |
|
321 |
xyzn3[0]=n3->get_x(); |
322 |
xyzn3[1]=n3->get_y(); |
323 |
xyzn3[2]=n3->get_z(); |
324 |
|
325 |
// Position des noeuds dans le repère initial (i) |
326 |
OT_VECTEUR_3D n1_i(xyzn1); |
327 |
OT_VECTEUR_3D n2_i(xyzn2); |
328 |
OT_VECTEUR_3D n3_i(xyzn3); |
329 |
OT_MATRICE_3D pos_i(n1_i,n2_i,n3_i); |
330 |
|
331 |
// Le nouveau repère est défini par les vecteurs vec1 et vec2 et son origine est le noeud n1 |
332 |
OT_VECTEUR_3D vec1(xyzn1,xyzn2); |
333 |
OT_VECTEUR_3D vec3(xyzn1,xyzn3); |
334 |
vec1.norme(); |
335 |
OT_VECTEUR_3D n=vec1&vec3; // Normale au triangle au noeud n1 |
336 |
n.norme(); |
337 |
OT_VECTEUR_3D vec2=n&vec1; // Vecteur normal à vec1 dans le plan du triangle |
338 |
vec2.norme(); |
339 |
OT_MATRICE_3D pif(vec1,vec2,n); // Matrice de passage du repère Ri au repère Rf |
340 |
|
341 |
// Position des noeuds dans le repère final (f) |
342 |
OT_MATRICE_3D pif_inverse=pif.inverse(); |
343 |
OT_MATRICE_3D pos_f=pif_inverse*pos_i; |
344 |
OT_VECTEUR_3D vecn1=pos_f.get_vecteur1(); |
345 |
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 |
346 |
pos_f.change_vecteur2(pos_f.get_vecteur2()-vecn1); |
347 |
pos_f.change_vecteur3(pos_f.get_vecteur3()-vecn1); |
348 |
|
349 |
// Déterminant du jacobien 2D dans le repère final |
350 |
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)); |
351 |
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)); |
352 |
OT_VECTEUR_3D dNi_dZ(0.,0.,0.); |
353 |
|
354 |
OT_MATRICE_3D dNi(dNi_de,dNi_dn,dNi_dZ); |
355 |
OT_MATRICE_3D dNi_transpose=dNi.transpose(); |
356 |
OT_MATRICE_3D pos_f_transpose=pos_f.transpose(); |
357 |
OT_MATRICE_3D J=dNi_transpose*pos_f_transpose; |
358 |
|
359 |
double J11=J(0,0); |
360 |
double J21=J(1,0); |
361 |
double J12=J(0,1); |
362 |
double J22=J(1,1); |
363 |
double det_J=J11*J22-J12*J21; |
364 |
|
365 |
// Jacobien inverse 3D dans le repère final |
366 |
OT_VECTEUR_3D col1_jf_3d(J22/det_J,-J21/det_J,0.); |
367 |
OT_VECTEUR_3D col2_jf_3d(-J12/det_J,J11/det_J,0.); |
368 |
OT_VECTEUR_3D col3_jf_3d(0.,0.,0.); |
369 |
OT_MATRICE_3D j_f_3d(col1_jf_3d,col2_jf_3d,col3_jf_3d); |
370 |
|
371 |
// Jacobien inverse 3D dans le repère initial |
372 |
OT_MATRICE_3D j_i=pif*j_f_3d; |
373 |
|
374 |
j[0]=j_i(0,0); |
375 |
j[1]=j_i(1,0); |
376 |
j[2]=j_i(2,0); |
377 |
|
378 |
j[3]=j_i(0,1); |
379 |
j[4]=j_i(1,1); |
380 |
j[5]=j_i(2,1); |
381 |
|
382 |
j[6]=j_i(0,2); |
383 |
j[7]=j_i(1,2); |
384 |
j[8]=j_i(2,2); |
385 |
} |