ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_triangle6.cpp
Revision: 786
Committed: Wed Mar 9 21:10:49 2016 UTC (9 years, 2 months ago) by francois
File size: 15602 byte(s)
Log Message: 
Possibilité d'obtenir la matrice de rigidite de code aster sous form de matrice bande.
Pour cela ajouter -opt dans l'operation -fem à la création du maillage FEM et dans -calculaster utiliser les nouveaux types d'étude rigidite.
Attention pour cela il faut une version de code aster modifiée
De manière automatique elle peut etre optenu en faisant ./instal_aster 12.5mod
Uniquement implantée avec la version 12.5 de code aster.

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_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 francois 612 #include "math.h"
31     #include "ot_quadrature_gauss.h"
32 francois 283
33 francois 309 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 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     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 francois 283 }
44    
45 francois 309 FEM_TRIANGLE6::FEM_TRIANGLE6(class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(mai),FEM_TEMPLATE_ELEMENT<6>(tabnoeud)
46 francois 283 {
47     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
48 francois 309 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 francois 283 }
56    
57 francois 378 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 francois 684 void FEM_TRIANGLE6::reinit_boite_3D(void)
107     {
108     FEM_TEMPLATE_ELEMENT<6>::reinit_boite_3D();
109     }
110 francois 378
111    
112    
113    
114 francois 309 FEM_TRIANGLE6::FEM_TRIANGLE6(FEM_TRIANGLE6& mdd):FEM_ELEMENT2(mdd),FEM_TEMPLATE_ELEMENT<6>(mdd)
115 francois 283 {
116     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
117 francois 309 tab[0]->get_lien_element2()->ajouter(this);
118     tab[1]->get_lien_element2()->ajouter(this);
119     tab[2]->get_lien_element2()->ajouter(this);
120     tab[3]->get_lien_element2()->ajouter(this);
121     tab[4]->get_lien_element2()->ajouter(this);
122     tab[5]->get_lien_element2()->ajouter(this);
123     get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
124 francois 283 }
125     FEM_TRIANGLE6::~FEM_TRIANGLE6()
126     {
127 francois 663 if (liaison_topologique!=NULL)
128     if (liaison_topologique->get_dimension()==0) liaison_topologique->get_lien_fem_maillage()->supprimer(this);
129 francois 309 tab[0]->get_lien_element2()->supprimer(this);
130     tab[1]->get_lien_element2()->supprimer(this);
131     tab[2]->get_lien_element2()->supprimer(this);
132     tab[3]->get_lien_element2()->supprimer(this);
133     tab[4]->get_lien_element2()->supprimer(this);
134     tab[5]->get_lien_element2()->supprimer(this);
135     get_fem_noeudpetitid()->get_lien_petit_element2()->supprimer(this);
136 francois 283 }
137    
138     FEM_ELEMENT_MAILLAGE* FEM_TRIANGLE6::dupliquer(FEM_MAILLAGE *femmai,long decalage)
139     {
140     FEM_NOEUD* tabnoeud[6];
141     tabnoeud[0]=femmai->get_fem_noeudid(tab[0]->get_id()+decalage);
142     tabnoeud[1]=femmai->get_fem_noeudid(tab[1]->get_id()+decalage);
143     tabnoeud[2]=femmai->get_fem_noeudid(tab[2]->get_id()+decalage);
144     tabnoeud[3]=femmai->get_fem_noeudid(tab[3]->get_id()+decalage);
145     tabnoeud[4]=femmai->get_fem_noeudid(tab[4]->get_id()+decalage);
146     tabnoeud[5]=femmai->get_fem_noeudid(tab[5]->get_id()+decalage);
147     FEM_TRIANGLE6* tri=new FEM_TRIANGLE6(get_id()+decalage,maillage,tabnoeud);
148 francois 309 femmai->ajouter_fem_element2(tri);
149 francois 283 return tri;
150     }
151    
152 francois 786
153     void FEM_TRIANGLE6::get_voisin_noeud(class FEM_NOEUD* no,TPL_LISTE_ENTITE<FEM_NOEUD*> &voisin)
154     {
155     voisin.vide();
156     if (no==tab[0])
157     {
158     voisin.ajouter(tab[1]);
159     voisin.ajouter(tab[5]);
160     }
161     if (no==tab[1])
162     {
163     voisin.ajouter(tab[0]);
164     voisin.ajouter(tab[2]);
165     }
166     if (no==tab[2])
167     {
168     voisin.ajouter(tab[1]);
169     voisin.ajouter(tab[3]);
170     }
171     if (no==tab[3])
172     {
173     voisin.ajouter(tab[2]);
174     voisin.ajouter(tab[4]);
175     }
176     if (no==tab[4])
177     {
178     voisin.ajouter(tab[3]);
179     voisin.ajouter(tab[5]);
180     }
181     if (no==tab[5])
182     {
183     voisin.ajouter(tab[4]);
184     voisin.ajouter(tab[0]);
185     }
186     }
187 francois 635 bool FEM_TRIANGLE6::valide_parametre_EF(double* uvw)
188     {
189     if (uvw[0]>=-1e-10)
190     if (uvw[1]>=-1e-10)
191     if (uvw[0]+uvw[1]<=1.+1e-10)
192     return true;
193     return false;
194     }
195 francois 283
196    
197     int FEM_TRIANGLE6::get_type_entite(void)
198     {
199     return IDFEM_TRIANGLE6;
200     }
201    
202     int FEM_TRIANGLE6::get_dimension(void)
203     {
204     return 2;
205     }
206    
207    
208     int FEM_TRIANGLE6::get_nb_fem_noeud(void)
209     {
210     return FEM_TEMPLATE_ELEMENT<6>::get_nb_fem_noeud();
211     }
212    
213     FEM_NOEUD* FEM_TRIANGLE6::get_fem_noeud(int num)
214     {
215     return FEM_TEMPLATE_ELEMENT<6>::get_fem_noeud(num);
216     }
217    
218     void FEM_TRIANGLE6::change_noeud(int num,FEM_NOEUD* noeud)
219     {
220     FEM_TEMPLATE_ELEMENT<6>::change_noeud(num,noeud);
221     }
222    
223     BOITE_3D& FEM_TRIANGLE6::get_boite_3D(void)
224     {
225     return FEM_TEMPLATE_ELEMENT<6>::get_boite_3D();
226     }
227    
228    
229 francois 763 void FEM_TRIANGLE6::enregistrer(std::ostream& o,double version)
230 francois 283 {
231 francois 378 if (maillage!=NULL)
232     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;
233     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;
234     else
235     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;
236     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;
237    
238 francois 283 }
239 francois 310
240 francois 637 int FEM_TRIANGLE6::get_nb_pt_gauss(int degre)
241     {
242 francois 754 return OT_POINTS_GAUSS::get_nb_point_tri(degre);
243 francois 637 }
244 francois 310
245 francois 637 void FEM_TRIANGLE6::get_pt_gauss(int degre,int num,double &w,double *uv)
246     {
247 francois 754 return OT_POINTS_GAUSS::get_pt_gauss_tri(degre,num,w,uv);
248 francois 637 }
249    
250 francois 757 int FEM_TRIANGLE6::get_degre_gauss(int num)
251     {
252     return OT_POINTS_GAUSS::get_degre_gauss_tri(num);
253     }
254    
255 francois 638 int FEM_TRIANGLE6::get_nb_fonction_interpolation(void)
256 francois 310 {
257     return 6;
258     }
259 francois 638 int FEM_TRIANGLE6::get_degremax_fonction_interpolation(void)
260     {
261     return 2;
262     }
263 francois 310 double FEM_TRIANGLE6::get_fonction_interpolation(int num,double *uv)
264     {
265     double val;
266     switch (num)
267     {
268     case 1:
269     val=(-1+uv[0]+uv[1])*(-1+2*uv[0]+2*uv[1]);
270     break;
271     case 2:
272     val=4.*uv[0]*(1-uv[0]-uv[1]);
273     break;
274     case 3:
275     val=-uv[0]*(1.-2.*uv[0]);
276     break;
277     case 4:
278     val=4.*uv[0]*uv[1];
279     break;
280     case 5:
281     val=-uv[1]*(1-2*uv[1]);
282     break;
283     case 6:
284     val=4.*uv[1]*(1-uv[0]-uv[1]);
285     break;
286    
287     }
288     return val;
289     }
290    
291     double FEM_TRIANGLE6::get_fonction_derive_interpolation(int num,int num_variable,double *uv)
292     {
293     double val;
294    
295     switch (num)
296     {
297     case 1:
298     switch (num_variable)
299     {
300     case 1:
301 francois 612 val=-3.+4*uv[0]+4*uv[1];
302 francois 310 break;
303     case 2:
304     val=-3.+4*uv[0]+4*uv[1];
305     break;
306     } break;
307     case 2:
308     switch (num_variable)
309     {
310     case 1:
311     val=4.-8*uv[0]-4*uv[1];
312     break;
313     case 2:
314     val=-4.*uv[0];
315     break;
316     }break;
317     case 3:
318     switch (num_variable)
319     {
320     case 1:
321 francois 612 val=-1.+4*uv[0];
322 francois 310 break;
323     case 2:
324     val=0.;
325     break;
326     }break;
327     case 4:
328     switch (num_variable)
329     {
330     case 1:
331     val=4*uv[1];
332     break;
333     case 2:
334     val=4*uv[0];
335     break;
336     }break;
337     case 5:
338     switch (num_variable)
339     {
340     case 1:
341     val=0.;
342     break;
343     case 2:
344     val=-1+4.*uv[1];
345     break;
346     }break;
347     case 6:
348     switch (num_variable)
349     {
350     case 1:
351     val=-4*uv[1];
352     break;
353     case 2:
354     val=4-4*uv[0]-8*uv[1];
355     break;
356     }break;
357    
358     }
359     return val;
360     }
361    
362 francois 612
363 francois 310 double FEM_TRIANGLE6::get_jacobien(double* jac,double *uv,int& li,int& col,double unite)
364     {
365 francois 612 FEM_NOEUD* n1=get_fem_noeud(0);
366     FEM_NOEUD* n2=get_fem_noeud(1);
367     FEM_NOEUD* n3=get_fem_noeud(2);
368     FEM_NOEUD* n4=get_fem_noeud(3);
369     FEM_NOEUD* n5=get_fem_noeud(4);
370     FEM_NOEUD* n6=get_fem_noeud(5);
371    
372     double xyzn1[3],xyzn2[3],xyzn3[3],xyzn4[3],xyzn5[3],xyzn6[3];
373    
374     xyzn1[0]=n1->get_x();
375     xyzn1[1]=n1->get_y();
376     xyzn1[2]=n1->get_z();
377    
378     xyzn2[0]=n2->get_x();
379     xyzn2[1]=n2->get_y();
380     xyzn2[2]=n2->get_z();
381    
382     xyzn3[0]=n3->get_x();
383     xyzn3[1]=n3->get_y();
384     xyzn3[2]=n3->get_z();
385    
386     xyzn4[0]=n4->get_x();
387     xyzn4[1]=n4->get_y();
388     xyzn4[2]=n4->get_z();
389    
390     xyzn5[0]=n5->get_x();
391     xyzn5[1]=n5->get_y();
392     xyzn5[2]=n5->get_z();
393    
394     xyzn6[0]=n6->get_x();
395     xyzn6[1]=n6->get_y();
396     xyzn6[2]=n6->get_z();
397    
398     OT_VECTEUR_3D vec1(xyzn1,xyzn3);
399     OT_VECTEUR_3D vec3(xyzn1,xyzn5);
400     OT_VECTEUR_3D n=vec1&vec3;
401     n.norme();
402 francois 310
403 francois 612 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];
404     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];
405     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];
406    
407     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];
408     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];
409     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];
410    
411     jac[6]=n[0];
412     jac[7]=n[1];
413     jac[8]=n[2];
414    
415     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]);
416    
417     return DET;
418    
419 francois 310 }
420    
421     void FEM_TRIANGLE6::get_inverse_jacob(double* j,double *uv,double unite)
422     {
423 francois 612 double jac[9];
424     int lig,col;
425     double detj=this->get_jacobien(jac,uv,lig,col,unite);
426    
427     j[0] = (jac[4]*jac[8]-jac[5]*jac[7])/(detj);
428     j[1] = -(jac[1]*jac[8]-jac[2]*jac[7])/(detj);
429     j[2] =-(-jac[1]*jac[5]+jac[2]*jac[4])/(detj);
430    
431     j[3] = -(jac[1]*jac[8]-jac[5]*jac[6])/(detj);
432     j[4] = (jac[0]*jac[8]-jac[2]*jac[6])/(detj);
433     j[5] = -(jac[0]*jac[5]-jac[2]*jac[3])/(detj);
434    
435     j[6] =-(-jac[3]*jac[7]+jac[4]*jac[6])/(detj);
436     j[7] = -(jac[0]*jac[7]-jac[1]*jac[6])/(detj);
437     j[8] = (jac[0]*jac[4]-jac[1]*jac[4])/(detj);
438    
439     }
440