ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_hexa8.cpp
Revision: 786
Committed: Wed Mar 9 21:10:49 2016 UTC (9 years, 2 months ago) by francois
File size: 16529 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 310 //------------------------------------------------------------
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_hexa8.h"
27     #include "fem_maillage.h"
28     #include "fem_noeud.h"
29     #include "mg_element_maillage.h"
30     #include "ot_tenseur.h"
31 francois 754 #include "ot_quadrature_gauss.h"
32 francois 310
33    
34     FEM_HEXA8::FEM_HEXA8(unsigned long num,class MG_ELEMENT_MAILLAGE* mai,class FEM_NOEUD** tabnoeud):FEM_ELEMENT3(num,mai),FEM_TEMPLATE_ELEMENT<8>(tabnoeud)
35     {
36 francois 581 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
37 francois 310 tab[0]->get_lien_element3()->ajouter(this);
38     tab[1]->get_lien_element3()->ajouter(this);
39     tab[2]->get_lien_element3()->ajouter(this);
40     tab[3]->get_lien_element3()->ajouter(this);
41     tab[4]->get_lien_element3()->ajouter(this);
42     tab[5]->get_lien_element3()->ajouter(this);
43     tab[6]->get_lien_element3()->ajouter(this);
44     tab[7]->get_lien_element3()->ajouter(this);
45     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
46     }
47    
48     FEM_HEXA8::FEM_HEXA8(class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT3(mai),FEM_TEMPLATE_ELEMENT<8>(tabnoeud)
49     {
50 francois 581 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
51 francois 310 tab[0]->get_lien_element3()->ajouter(this);
52     tab[1]->get_lien_element3()->ajouter(this);
53     tab[2]->get_lien_element3()->ajouter(this);
54     tab[3]->get_lien_element3()->ajouter(this);
55     tab[4]->get_lien_element3()->ajouter(this);
56     tab[5]->get_lien_element3()->ajouter(this);
57     tab[6]->get_lien_element3()->ajouter(this);
58     tab[7]->get_lien_element3()->ajouter(this);
59     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
60     }
61    
62 francois 378 FEM_HEXA8::FEM_HEXA8(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo,class FEM_NOEUD** tabnoeud):FEM_ELEMENT3(num,topo),FEM_TEMPLATE_ELEMENT<8>(tabnoeud)
63     {
64 francois 581 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
65 francois 378 tab[0]->get_lien_element3()->ajouter(this);
66     tab[1]->get_lien_element3()->ajouter(this);
67     tab[2]->get_lien_element3()->ajouter(this);
68     tab[3]->get_lien_element3()->ajouter(this);
69     tab[4]->get_lien_element3()->ajouter(this);
70     tab[5]->get_lien_element3()->ajouter(this);
71     tab[6]->get_lien_element3()->ajouter(this);
72     tab[7]->get_lien_element3()->ajouter(this);
73     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
74     }
75    
76     FEM_HEXA8::FEM_HEXA8(class MG_ELEMENT_TOPOLOGIQUE* topo,FEM_NOEUD** tabnoeud):FEM_ELEMENT3(topo),FEM_TEMPLATE_ELEMENT<8>(tabnoeud)
77     {
78 francois 581 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
79 francois 378 tab[0]->get_lien_element3()->ajouter(this);
80     tab[1]->get_lien_element3()->ajouter(this);
81     tab[2]->get_lien_element3()->ajouter(this);
82     tab[3]->get_lien_element3()->ajouter(this);
83     tab[4]->get_lien_element3()->ajouter(this);
84     tab[5]->get_lien_element3()->ajouter(this);
85     tab[6]->get_lien_element3()->ajouter(this);
86     tab[7]->get_lien_element3()->ajouter(this);
87     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
88     }
89    
90     FEM_HEXA8::FEM_HEXA8(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai,class FEM_NOEUD** tabnoeud):FEM_ELEMENT3(num,topo,mai),FEM_TEMPLATE_ELEMENT<8>(tabnoeud)
91     {
92 francois 581 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
93 francois 378 tab[0]->get_lien_element3()->ajouter(this);
94     tab[1]->get_lien_element3()->ajouter(this);
95     tab[2]->get_lien_element3()->ajouter(this);
96     tab[3]->get_lien_element3()->ajouter(this);
97     tab[4]->get_lien_element3()->ajouter(this);
98     tab[5]->get_lien_element3()->ajouter(this);
99     tab[6]->get_lien_element3()->ajouter(this);
100     tab[7]->get_lien_element3()->ajouter(this);
101     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
102     }
103    
104     FEM_HEXA8::FEM_HEXA8(class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT3(topo,mai),FEM_TEMPLATE_ELEMENT<8>(tabnoeud)
105     {
106 francois 581 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
107 francois 378 tab[0]->get_lien_element3()->ajouter(this);
108     tab[1]->get_lien_element3()->ajouter(this);
109     tab[2]->get_lien_element3()->ajouter(this);
110     tab[3]->get_lien_element3()->ajouter(this);
111     tab[4]->get_lien_element3()->ajouter(this);
112     tab[5]->get_lien_element3()->ajouter(this);
113     tab[6]->get_lien_element3()->ajouter(this);
114     tab[7]->get_lien_element3()->ajouter(this);
115     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
116     }
117    
118    
119    
120    
121 francois 310 FEM_HEXA8::FEM_HEXA8(FEM_HEXA8& mdd):FEM_ELEMENT3(mdd),FEM_TEMPLATE_ELEMENT<8>(mdd)
122     {
123 francois 581 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
124 francois 310 tab[0]->get_lien_element3()->ajouter(this);
125     tab[1]->get_lien_element3()->ajouter(this);
126     tab[2]->get_lien_element3()->ajouter(this);
127     tab[3]->get_lien_element3()->ajouter(this);
128     tab[4]->get_lien_element3()->ajouter(this);
129     tab[5]->get_lien_element3()->ajouter(this);
130     tab[6]->get_lien_element3()->ajouter(this);
131     tab[7]->get_lien_element3()->ajouter(this);
132     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
133     }
134     FEM_HEXA8::~FEM_HEXA8()
135     {
136 francois 663 if (liaison_topologique!=NULL)
137     if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->supprimer(this);
138 francois 310 tab[0]->get_lien_element3()->supprimer(this);
139     tab[1]->get_lien_element3()->supprimer(this);
140     tab[2]->get_lien_element3()->supprimer(this);
141     tab[3]->get_lien_element3()->supprimer(this);
142     tab[4]->get_lien_element3()->supprimer(this);
143     tab[5]->get_lien_element3()->supprimer(this);
144     tab[6]->get_lien_element3()->supprimer(this);
145     tab[7]->get_lien_element3()->supprimer(this);
146     get_fem_noeudpetitid()->get_lien_petit_element3()->supprimer(this);
147     }
148    
149    
150     FEM_ELEMENT_MAILLAGE* FEM_HEXA8::dupliquer(FEM_MAILLAGE *femmai,long decalage)
151     {
152     FEM_NOEUD* tabnoeud[8];
153     tabnoeud[0]=femmai->get_fem_noeudid(tab[0]->get_id()+decalage);
154     tabnoeud[1]=femmai->get_fem_noeudid(tab[1]->get_id()+decalage);
155     tabnoeud[2]=femmai->get_fem_noeudid(tab[2]->get_id()+decalage);
156     tabnoeud[3]=femmai->get_fem_noeudid(tab[3]->get_id()+decalage);
157     tabnoeud[4]=femmai->get_fem_noeudid(tab[4]->get_id()+decalage);
158     tabnoeud[5]=femmai->get_fem_noeudid(tab[5]->get_id()+decalage);
159     tabnoeud[6]=femmai->get_fem_noeudid(tab[6]->get_id()+decalage);
160     tabnoeud[7]=femmai->get_fem_noeudid(tab[7]->get_id()+decalage);
161     FEM_HEXA8* quad=new FEM_HEXA8(get_id()+decalage,maillage,tabnoeud);
162     femmai->ajouter_fem_element3(quad);
163     return quad;
164     }
165 francois 786 void FEM_HEXA8::get_voisin_noeud(class FEM_NOEUD* no,TPL_LISTE_ENTITE<FEM_NOEUD*> &voisin)
166     {
167     voisin.vide();
168     if (no==tab[0])
169     {
170     voisin.ajouter(tab[1]);
171     voisin.ajouter(tab[3]);
172     voisin.ajouter(tab[4]);
173     }
174     if (no==tab[1])
175     {
176     voisin.ajouter(tab[0]);
177     voisin.ajouter(tab[2]);
178     voisin.ajouter(tab[5]);
179     }
180     if (no==tab[2])
181     {
182     voisin.ajouter(tab[1]);
183     voisin.ajouter(tab[3]);
184     voisin.ajouter(tab[6]);
185     }
186     if (no==tab[3])
187     {
188     voisin.ajouter(tab[0]);
189     voisin.ajouter(tab[2]);
190     voisin.ajouter(tab[7]);
191     }
192     if (no==tab[4])
193     {
194     voisin.ajouter(tab[0]);
195     voisin.ajouter(tab[5]);
196     voisin.ajouter(tab[7]);
197     }
198     if (no==tab[5])
199     {
200     voisin.ajouter(tab[1]);
201     voisin.ajouter(tab[4]);
202     voisin.ajouter(tab[6]);
203     }
204     if (no==tab[6])
205     {
206     voisin.ajouter(tab[2]);
207     voisin.ajouter(tab[5]);
208     voisin.ajouter(tab[7]);
209     }
210     if (no==tab[7])
211     {
212     voisin.ajouter(tab[3]);
213     voisin.ajouter(tab[4]);
214     voisin.ajouter(tab[6]);
215     }
216     }
217 francois 310
218 francois 635 bool FEM_HEXA8::valide_parametre_EF(double* uvw)
219     {
220 francois 674 if (uvw[0]>=-1.-1e-10)
221     if (uvw[1]>=-1.-1e-10)
222     if (uvw[2]>=-1.-1e-10)
223 francois 635 if (uvw[0]<=1.+1e-10)
224     if (uvw[1]<=1.+1e-10)
225     if (uvw[2]<=1.+1e-10)
226     return true;
227     return false;
228     }
229 francois 310
230 francois 684 void FEM_HEXA8::reinit_boite_3D(void)
231     {
232     FEM_TEMPLATE_ELEMENT<8>::reinit_boite_3D();
233     }
234 francois 310
235     int FEM_HEXA8::get_type_entite(void)
236     {
237     return IDFEM_HEXA8;
238     }
239    
240     int FEM_HEXA8::get_dimension(void)
241     {
242     return 3;
243     }
244    
245    
246    
247     int FEM_HEXA8::get_nb_fem_noeud(void)
248     {
249     return FEM_TEMPLATE_ELEMENT<8>::get_nb_fem_noeud();
250     }
251    
252     FEM_NOEUD* FEM_HEXA8::get_fem_noeud(int num)
253     {
254     return FEM_TEMPLATE_ELEMENT<8>::get_fem_noeud(num);
255     }
256    
257     void FEM_HEXA8::change_noeud(int num,FEM_NOEUD* noeud)
258     {
259     FEM_TEMPLATE_ELEMENT<8>::change_noeud(num,noeud);
260     }
261    
262     BOITE_3D& FEM_HEXA8::get_boite_3D(void)
263     {
264     return FEM_TEMPLATE_ELEMENT<8>::get_boite_3D();
265     }
266    
267    
268 francois 763 void FEM_HEXA8::enregistrer(std::ostream& o,double version)
269 francois 310 {
270 francois 378 if (maillage!=NULL)
271     if (get_lien_topologie()!=NULL) o << "%" << get_id() << "=FEM_HEXA8($"<< 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() << ",$" << tab[6]->get_id()<< ",$" << tab[7]->get_id()<< ");" << std::endl;
272     else o << "%" << get_id() << "=FEM_HEXA8(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() << ",$" << tab[6]->get_id()<< ",$" << tab[7]->get_id()<<");" << std::endl;
273     else
274     if (get_lien_topologie()!=NULL) o << "%" << get_id() << "=FEM_HEXA8($"<< 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() << ",$" << tab[6]->get_id()<< ",$" << tab[7]->get_id()<< ");" << std::endl;
275     else o << "%" << get_id() << "=FEM_HEXA8(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() << ",$" << tab[6]->get_id()<< ",$" << tab[7]->get_id()<<");" << std::endl;
276    
277 francois 310 }
278 francois 637
279    
280     int FEM_HEXA8::get_nb_pt_gauss(int degre)
281     {
282 francois 754 return OT_POINTS_GAUSS::get_nb_point_hexa_prod(degre);
283 francois 637 }
284    
285     void FEM_HEXA8::get_pt_gauss(int degre,int num,double &w,double *uvw)
286     {
287 francois 754 return OT_POINTS_GAUSS::get_pt_gauss_hex_prod(degre,num,w,uvw);
288 francois 637 }
289 francois 754
290 francois 757 int FEM_HEXA8::get_degre_gauss(int num)
291     {
292     return OT_POINTS_GAUSS::get_degre_gauss_hexa_prod(num);
293     }
294    
295 francois 638 int FEM_HEXA8::get_nb_fonction_interpolation(void)
296 francois 310 {
297     return 8;
298     }
299    
300 francois 638 int FEM_HEXA8::get_degremax_fonction_interpolation(void)
301     {
302     return 3;
303     }
304 francois 310 double FEM_HEXA8::get_fonction_interpolation(int num,double *uv)
305     {
306     double val;
307     double a1=1.+uv[0],a2=1.-uv[0];
308     double b1=1.+uv[1],b2=1.-uv[1];
309     double c1=1.+uv[2],c2=1.-uv[2];
310     switch (num)
311     {
312     case 1:
313     val=0.125*a2*b2*c2;
314     break;
315     case 2:
316     val=0.125*a1*b2*c2;
317     break;
318     case 3:
319     val=0.125*a1*b1*c2;
320     break;
321     case 4:
322     val=0.125*a2*b1*c2;
323     break;
324     case 5:
325     val=0.125*a2*b2*c1;
326     break;
327     case 6:
328     val=0.125*a1*b2*c1;
329     break;
330     case 7:
331     val=0.125*a1*b1*c1;
332     break;
333     case 8:
334     val=0.125*a2*b1*c1;
335     break;
336     }
337     return val;
338     }
339    
340     double FEM_HEXA8::get_fonction_derive_interpolation(int num,int num_variable,double *uv)
341     {
342     double val;
343     double a1=1.+uv[0],a2=1.-uv[0];
344     double b1=1.+uv[1],b2=1.-uv[1];
345     double c1=1.+uv[2],c2=1.-uv[2];
346     switch (num)
347     {
348     case 1:
349     switch (num_variable)
350     {
351     case 1:
352     val=-0.125*b2*c2;
353     break;
354     case 2:
355     val=-0.125*a2*c2;
356     break;
357     case 3:
358     val=-0.125*a2*b2;
359     break;
360     } break;
361     case 2:
362     switch (num_variable)
363     {
364     case 1:
365     val=0.125*b2*c2;
366     break;
367     case 2:
368     val=-0.125*a1*c2;
369     break;
370     case 3:
371     val=-0.125*a1*b2;
372     break;
373     } break;
374     case 3:
375     switch (num_variable)
376     {
377     case 1:
378     val=0.125*b1*c2;
379     break;
380     case 2:
381     val=0.125*a1*c2;
382     break;
383     case 3:
384     val=-0.125*a1*b1;
385     break;
386     } break;
387     case 4:
388     switch (num_variable)
389     {
390     case 1:
391     val=-0.125*b1*c2;
392     break;
393     case 2:
394     val=0.125*a2*c2;
395     break;
396     case 3:
397     val=-0.125*a2*b1;
398     break;
399     } break;
400     case 5:
401     switch (num_variable)
402     {
403     case 1:
404     val=-0.125*b2*c1;
405     break;
406     case 2:
407     val=-0.125*a2*c1;
408     break;
409     case 3:
410     val=0.125*a2*b2;
411     break;
412     } break;
413     case 6:
414     switch (num_variable)
415     {
416     case 1:
417     val=0.125*b2*c1;
418     break;
419     case 2:
420     val=-0.125*a1*c1;
421     break;
422     case 3:
423     val=0.125*a1*b2;
424     break;
425     } break;
426     case 7:
427     switch (num_variable)
428     {
429     case 1:
430     val=0.125*b1*c1;
431     break;
432     case 2:
433     val=0.125*a1*c1;
434     break;
435     case 3:
436     val=0.125*a1*b1;
437     break;
438     } break;
439     case 8:
440     switch (num_variable)
441     {
442     case 1:
443     val=-0.125*b1*c1;
444     break;
445     case 2:
446     val=0.125*a2*c1;
447     break;
448     case 3:
449     val=0.125*a2*b1;
450     break;
451     } break;
452     }
453     return val;
454     }
455    
456     double FEM_HEXA8::get_jacobien(double* jac,double *uv,int& li,int& col,double unite)
457     {
458     li=3;col=3;
459     OT_TENSEUR Ni(3,8),Xi(8,3);
460     for (int i=0;i<3;i++)
461     for (int j=0;j<8;j++)
462     Ni(i,j)=get_fonction_derive_interpolation(j+1,i+1,uv);
463     for (int i=0;i<8;i++)
464     {
465     double *xyz=tab[i]->get_coord();
466     Xi(i,0)=xyz[0]*unite;
467     Xi(i,1)=xyz[1]*unite;
468     Xi(i,2)=xyz[2]*unite;
469     }
470     OT_TENSEUR jacobi=Ni*Xi;
471     jac[0]=jacobi(0,0).get_x();
472     jac[1]=jacobi(0,1).get_x();
473     jac[2]=jacobi(0,2).get_x();
474     jac[3]=jacobi(1,0).get_x();
475     jac[4]=jacobi(1,1).get_x();
476     jac[5]=jacobi(1,2).get_x();
477     jac[6]=jacobi(2,0).get_x();
478     jac[7]=jacobi(2,1).get_x();
479     jac[8]=jacobi(2,2).get_x();
480     double det=jac[0]*jac[4]*jac[8]+jac[3]*jac[7]*jac[2]+jac[6]*jac[1]*jac[5]-jac[2]*jac[4]*jac[6]-jac[5]*jac[7]*jac[0]-jac[8]*jac[1]*jac[3];
481     return det;
482    
483     }
484    
485     void FEM_HEXA8::get_inverse_jacob(double* j,double *uv,double unite)
486     {
487     double jac[9];
488     int li,col;
489     double detj=get_jacobien(jac,uv,li,col,unite);
490     j[0*3+0] = (jac[1*3+1]*jac[2*3+2]-jac[1*3+2]*jac[2*3+1])/(detj);
491     j[0*3+1] = -(jac[0*3+1]*jac[2*3+2]-jac[0*3+2]*jac[2*3+1])/(detj);
492     j[0*3+2] =-(-jac[0*3+1]*jac[1*3+2]+jac[0*3+2]*jac[1*3+1])/(detj);
493    
494     j[1*3+0] = -(jac[1*3+0]*jac[2*3+2]-jac[1*3+2]*jac[2*3+0])/(detj);
495     j[1*3+1] = (jac[0*3+0]*jac[2*3+2]-jac[0*3+2]*jac[2*3+0])/(detj);
496     j[1*3+2] = -(jac[0*3+0]*jac[1*3+2]-jac[0*3+2]*jac[1*3+0])/(detj);
497    
498     j[2*3+0] =-(-jac[1*3+0]*jac[2*3+1]+jac[1*3+1]*jac[2*3+0])/(detj);
499     j[2*3+1] = -(jac[0*3+0]*jac[2*3+1]-jac[0*3+1]*jac[2*3+0])/(detj);
500     j[2*3+2] = (jac[0*3+0]*jac[1*3+1]-jac[0*3+1]*jac[1*3+0])/(detj);
501    
502 francois 405 }
503    
504    
505    
506 francois 406 int FEM_HEXA8::verifie_validite_decoupage_xfem(double *vol)
507 francois 405 {
508    
509 francois 310 }