ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_triangle3.cpp
Revision: 392
Committed: Mon Mar 11 20:35:58 2013 UTC (12 years, 5 months ago) by gervaislavoie
File size: 14055 byte(s)
Log Message:
Corrections calcul jacobien/jacobien inverse et ajout jacobien inverse dans opt_triangle.

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::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 }