ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_triangle3.cpp
Revision: 389
Committed: Tue Feb 12 17:06:17 2013 UTC (12 years, 3 months ago) by gervaislavoie
File size: 14321 byte(s)
Log Message: 
Probleme lors de la mise a jour automatique pour le calcul de l'angle, donc fait manuellement.
Calcul du jacobien et jacobien inverse a partir des coordonnees initiales des noeuds.
Classification du code en fonctions pour en faciliter la comprehension.

File Contents

# User Rev Content
1 francois 283 //------------------------------------------------------------
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 gervaislavoie 382 #include "math.h"
31 francois 283
32    
33 francois 309 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 francois 283 {
35     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
36 francois 309 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 francois 283 }
41    
42 francois 309 FEM_TRIANGLE3::FEM_TRIANGLE3(class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(mai),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
43 francois 283 {
44     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
45 francois 309 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 francois 283 }
50 francois 378 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 francois 283
59 francois 378 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 francois 309 FEM_TRIANGLE3::FEM_TRIANGLE3(FEM_TRIANGLE3& mdd):FEM_ELEMENT2(mdd),FEM_TEMPLATE_ELEMENT<3>(mdd)
85 francois 283 {
86     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
87 francois 309 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 francois 283 }
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 francois 309 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 francois 283 }
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 francois 309 femmai->ajouter_fem_element2(tri);
111 francois 283 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 francois 378 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 francois 310
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 gervaislavoie 382 {
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 gervaislavoie 389 // 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 gervaislavoie 382 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 gervaislavoie 389 double xyzn1[3];
236     double xyzn2[3];
237     double xyzn3[3];
238    
239     xyzn1[0]=n1->get_dx(); // Position initiale
240     xyzn1[1]=n1->get_dy();
241     xyzn1[2]=n1->get_dz();
242    
243     xyzn2[0]=n2->get_dx();
244     xyzn2[1]=n2->get_dy();
245     xyzn2[2]=n2->get_dz();
246    
247     xyzn3[0]=n3->get_dx();
248     xyzn3[1]=n3->get_dy();
249     xyzn3[2]=n3->get_dz();
250    
251 gervaislavoie 382 // Le nouveau repère est défini par les vecteurs vec1 et vec2 et son origine est le noeud n1
252     OT_VECTEUR_3D vec1(xyzn1,xyzn2);
253     OT_VECTEUR_3D vec3(xyzn1,xyzn3);
254     vec1.norme();
255     OT_VECTEUR_3D n=vec1&vec3; // Normale au triangle au noeud n1
256     n.norme();
257     OT_VECTEUR_3D vec2=n&vec1; // Vecteur normal à vec1 dans le plan du triangle
258     vec2.norme();
259     OT_MATRICE_3D pif(vec1,vec2,n); // Matrice de passage du repère Ri au repère Rf
260    
261     // Position initiale (0) des noeuds dans le repère initial (i)
262     OT_VECTEUR_3D n1_0i(n1->get_dx(),n1->get_dy(),n1->get_dz());
263     OT_VECTEUR_3D n2_0i(n2->get_dx(),n2->get_dy(),n2->get_dz());
264     OT_VECTEUR_3D n3_0i(n3->get_dx(),n3->get_dy(),n3->get_dz());
265     OT_MATRICE_3D pos_0i(n1_0i,n2_0i,n3_0i);
266    
267     // Position initiale (0) des noeuds dans le repère final (f)
268     OT_MATRICE_3D pif_inverse=pif.inverse();
269     OT_MATRICE_3D pos_0f=pif_inverse*pos_0i;
270     OT_VECTEUR_3D vecn1_0=pos_0f.get_vecteur1();
271     pos_0f.change_vecteur1(pos_0f.get_vecteur1()-vecn1_0); // Étant donné que le noeud 1 est l'origine du repère final, on soustrait ses coordonnées à chaque noeud
272     pos_0f.change_vecteur2(pos_0f.get_vecteur2()-vecn1_0);
273     pos_0f.change_vecteur3(pos_0f.get_vecteur3()-vecn1_0);
274    
275 gervaislavoie 389 // 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 gervaislavoie 382
280 gervaislavoie 389 OT_MATRICE_3D dNi(dNi_de,dNi_dn,dNi_dZ);
281     OT_MATRICE_3D dNi_transpose=dNi.transpose();
282     OT_MATRICE_3D pos_0f_transpose=pos_0f.transpose();
283     OT_MATRICE_3D J_f=dNi_transpose*pos_0f_transpose;
284 gervaislavoie 382
285 gervaislavoie 389 // Jacobien dans le repère initial
286     OT_MATRICE_3D J_i=pif*J_f;
287 gervaislavoie 382
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 francois 310 }
300    
301     void FEM_TRIANGLE3::get_inverse_jacob(double* j,double *uv,double unite)
302     {
303 gervaislavoie 382 // 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 gervaislavoie 389 // x',y',z' où z' est nul afin de pouvoir utiliser les fonctions d'interpolation de l'élément de référence
305 gervaislavoie 382 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 gervaislavoie 389 double xyzn1[3];
310     double xyzn2[3];
311     double xyzn3[3];
312 gervaislavoie 382
313 gervaislavoie 389 xyzn1[0]=n1->get_dx(); // Position initiale
314     xyzn1[1]=n1->get_dy();
315     xyzn1[2]=n1->get_dz();
316 gervaislavoie 382
317 gervaislavoie 389 xyzn2[0]=n2->get_dx();
318     xyzn2[1]=n2->get_dy();
319     xyzn2[2]=n2->get_dz();
320 gervaislavoie 382
321 gervaislavoie 389 xyzn3[0]=n3->get_dx();
322     xyzn3[1]=n3->get_dy();
323     xyzn3[2]=n3->get_dz();
324 gervaislavoie 382
325     OT_VECTEUR_3D vec1(xyzn1,xyzn2);
326     OT_VECTEUR_3D vec3(xyzn1,xyzn3);
327     vec1.norme();
328     OT_VECTEUR_3D n=vec1&vec3; // Normale au triangle au noeud n1
329     n.norme();
330     OT_VECTEUR_3D vec2=n&vec1; // Vecteur normal à vec1 dans le plan du triangle
331     vec2.norme();
332     OT_MATRICE_3D pif(vec1,vec2,n); // Matrice de passage du repère Ri au repère Rf
333    
334     // Position initiale (0) des noeuds dans le repère initial (i)
335     OT_VECTEUR_3D n1_0i(n1->get_dx(),n1->get_dy(),n1->get_dz());
336     OT_VECTEUR_3D n2_0i(n2->get_dx(),n2->get_dy(),n2->get_dz());
337     OT_VECTEUR_3D n3_0i(n3->get_dx(),n3->get_dy(),n3->get_dz());
338     OT_MATRICE_3D pos_0i(n1_0i,n2_0i,n3_0i);
339    
340     // Position initiale (0) des noeuds dans le repère final (f)
341     OT_MATRICE_3D pif_inverse=pif.inverse();
342     OT_MATRICE_3D pos_0f=pif_inverse*pos_0i;
343     OT_VECTEUR_3D vecn1_0=pos_0f.get_vecteur1();
344     pos_0f.change_vecteur1(pos_0f.get_vecteur1()-vecn1_0); // Étant donné que le noeud 1 est l'origine du repère final, on soustrait ses coordonnées à chaque noeud
345     pos_0f.change_vecteur2(pos_0f.get_vecteur2()-vecn1_0);
346     pos_0f.change_vecteur3(pos_0f.get_vecteur3()-vecn1_0);
347    
348 gervaislavoie 389 // Déterminant du jacobien 2D dans le repère final
349     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));
350     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));
351     OT_VECTEUR_3D dNi_dZ(0.,0.,0.);
352 gervaislavoie 382
353 gervaislavoie 389 OT_MATRICE_3D dNi(dNi_de,dNi_dn,dNi_dZ);
354     OT_MATRICE_3D dNi_transpose=dNi.transpose();
355     OT_MATRICE_3D pos_0f_transpose=pos_0f.transpose();
356     OT_MATRICE_3D J=dNi_transpose*pos_0f_transpose;
357 gervaislavoie 382
358 gervaislavoie 389 double J11=J(0,0);
359     double J21=J(1,0);
360     double J12=J(0,1);
361     double J22=J(1,1);
362 gervaislavoie 382 double det_J=J11*J22-J12*J21;
363    
364 gervaislavoie 389 // Jacobien inverse 3D dans le repère final
365 gervaislavoie 382 OT_VECTEUR_3D col1_jf_3d(J22/det_J,-J21/det_J,0.);
366     OT_VECTEUR_3D col2_jf_3d(-J12/det_J,J11/det_J,0.);
367     OT_VECTEUR_3D col3_jf_3d(0.,0.,0.);
368 gervaislavoie 389 OT_MATRICE_3D j_f_3d(col1_jf_3d,col2_jf_3d,col3_jf_3d);
369 gervaislavoie 382
370 gervaislavoie 389 // Jacobien inverse 3D dans le repère initial
371     OT_MATRICE_3D j_i=pif*j_f_3d;
372 gervaislavoie 382
373     j[0]=j_i(0,0);
374     j[1]=j_i(1,0);
375     j[2]=j_i(2,0);
376    
377     j[3]=j_i(0,1);
378     j[4]=j_i(1,1);
379     j[5]=j_i(2,1);
380    
381     j[6]=j_i(0,2);
382     j[7]=j_i(1,2);
383     j[8]=j_i(2,2);
384 francois 310 }