ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_hexa8.cpp
Revision: 310
Committed: Thu Jan 26 21:14:45 2012 UTC (13 years, 3 months ago) by francois
File size: 10789 byte(s)
Log Message:
Ajout des elements quadrangle et hexaedre dans les maillages FEM
Changement de la version des fichier de visu sous GMSH. Les .pos ont disparus et tout est intégré dans les .msh

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    
32    
33     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)
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_element3()->ajouter(this);
37     tab[1]->get_lien_element3()->ajouter(this);
38     tab[2]->get_lien_element3()->ajouter(this);
39     tab[3]->get_lien_element3()->ajouter(this);
40     tab[4]->get_lien_element3()->ajouter(this);
41     tab[5]->get_lien_element3()->ajouter(this);
42     tab[6]->get_lien_element3()->ajouter(this);
43     tab[7]->get_lien_element3()->ajouter(this);
44     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
45     }
46    
47     FEM_HEXA8::FEM_HEXA8(class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT3(mai),FEM_TEMPLATE_ELEMENT<8>(tabnoeud)
48     {
49     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
50     tab[0]->get_lien_element3()->ajouter(this);
51     tab[1]->get_lien_element3()->ajouter(this);
52     tab[2]->get_lien_element3()->ajouter(this);
53     tab[3]->get_lien_element3()->ajouter(this);
54     tab[4]->get_lien_element3()->ajouter(this);
55     tab[5]->get_lien_element3()->ajouter(this);
56     tab[6]->get_lien_element3()->ajouter(this);
57     tab[7]->get_lien_element3()->ajouter(this);
58     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
59     }
60    
61     FEM_HEXA8::FEM_HEXA8(FEM_HEXA8& mdd):FEM_ELEMENT3(mdd),FEM_TEMPLATE_ELEMENT<8>(mdd)
62     {
63     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
64     tab[0]->get_lien_element3()->ajouter(this);
65     tab[1]->get_lien_element3()->ajouter(this);
66     tab[2]->get_lien_element3()->ajouter(this);
67     tab[3]->get_lien_element3()->ajouter(this);
68     tab[4]->get_lien_element3()->ajouter(this);
69     tab[5]->get_lien_element3()->ajouter(this);
70     tab[6]->get_lien_element3()->ajouter(this);
71     tab[7]->get_lien_element3()->ajouter(this);
72     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
73     }
74     FEM_HEXA8::~FEM_HEXA8()
75     {
76     if (liaison_topologique==NULL) return;
77     if (liaison_topologique->get_dimension()==0) liaison_topologique->get_lien_fem_maillage()->supprimer(this);
78     tab[0]->get_lien_element3()->supprimer(this);
79     tab[1]->get_lien_element3()->supprimer(this);
80     tab[2]->get_lien_element3()->supprimer(this);
81     tab[3]->get_lien_element3()->supprimer(this);
82     tab[4]->get_lien_element3()->supprimer(this);
83     tab[5]->get_lien_element3()->supprimer(this);
84     tab[6]->get_lien_element3()->supprimer(this);
85     tab[7]->get_lien_element3()->supprimer(this);
86     get_fem_noeudpetitid()->get_lien_petit_element3()->supprimer(this);
87     }
88    
89    
90     FEM_ELEMENT_MAILLAGE* FEM_HEXA8::dupliquer(FEM_MAILLAGE *femmai,long decalage)
91     {
92     FEM_NOEUD* tabnoeud[8];
93     tabnoeud[0]=femmai->get_fem_noeudid(tab[0]->get_id()+decalage);
94     tabnoeud[1]=femmai->get_fem_noeudid(tab[1]->get_id()+decalage);
95     tabnoeud[2]=femmai->get_fem_noeudid(tab[2]->get_id()+decalage);
96     tabnoeud[3]=femmai->get_fem_noeudid(tab[3]->get_id()+decalage);
97     tabnoeud[4]=femmai->get_fem_noeudid(tab[4]->get_id()+decalage);
98     tabnoeud[5]=femmai->get_fem_noeudid(tab[5]->get_id()+decalage);
99     tabnoeud[6]=femmai->get_fem_noeudid(tab[6]->get_id()+decalage);
100     tabnoeud[7]=femmai->get_fem_noeudid(tab[7]->get_id()+decalage);
101     FEM_HEXA8* quad=new FEM_HEXA8(get_id()+decalage,maillage,tabnoeud);
102     femmai->ajouter_fem_element3(quad);
103     return quad;
104     }
105    
106    
107    
108    
109    
110     int FEM_HEXA8::get_type_entite(void)
111     {
112     return IDFEM_HEXA8;
113     }
114    
115     int FEM_HEXA8::get_dimension(void)
116     {
117     return 3;
118     }
119    
120    
121    
122     int FEM_HEXA8::get_nb_fem_noeud(void)
123     {
124     return FEM_TEMPLATE_ELEMENT<8>::get_nb_fem_noeud();
125     }
126    
127     FEM_NOEUD* FEM_HEXA8::get_fem_noeud(int num)
128     {
129     return FEM_TEMPLATE_ELEMENT<8>::get_fem_noeud(num);
130     }
131    
132     void FEM_HEXA8::change_noeud(int num,FEM_NOEUD* noeud)
133     {
134     FEM_TEMPLATE_ELEMENT<8>::change_noeud(num,noeud);
135     }
136    
137     BOITE_3D& FEM_HEXA8::get_boite_3D(void)
138     {
139     return FEM_TEMPLATE_ELEMENT<8>::get_boite_3D();
140     }
141    
142    
143     void FEM_HEXA8::enregistrer(std::ostream& o)
144     {
145     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;
146     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;
147     }
148     int FEM_HEXA8::nb_fonction_interpolation(void)
149     {
150     return 8;
151     }
152    
153    
154     double FEM_HEXA8::get_fonction_interpolation(int num,double *uv)
155     {
156     double val;
157     double a1=1.+uv[0],a2=1.-uv[0];
158     double b1=1.+uv[1],b2=1.-uv[1];
159     double c1=1.+uv[2],c2=1.-uv[2];
160     switch (num)
161     {
162     case 1:
163     val=0.125*a2*b2*c2;
164     break;
165     case 2:
166     val=0.125*a1*b2*c2;
167     break;
168     case 3:
169     val=0.125*a1*b1*c2;
170     break;
171     case 4:
172     val=0.125*a2*b1*c2;
173     break;
174     case 5:
175     val=0.125*a2*b2*c1;
176     break;
177     case 6:
178     val=0.125*a1*b2*c1;
179     break;
180     case 7:
181     val=0.125*a1*b1*c1;
182     break;
183     case 8:
184     val=0.125*a2*b1*c1;
185     break;
186     }
187     return val;
188     }
189    
190     double FEM_HEXA8::get_fonction_derive_interpolation(int num,int num_variable,double *uv)
191     {
192     double val;
193     double a1=1.+uv[0],a2=1.-uv[0];
194     double b1=1.+uv[1],b2=1.-uv[1];
195     double c1=1.+uv[2],c2=1.-uv[2];
196     switch (num)
197     {
198     case 1:
199     switch (num_variable)
200     {
201     case 1:
202     val=-0.125*b2*c2;
203     break;
204     case 2:
205     val=-0.125*a2*c2;
206     break;
207     case 3:
208     val=-0.125*a2*b2;
209     break;
210     } break;
211     case 2:
212     switch (num_variable)
213     {
214     case 1:
215     val=0.125*b2*c2;
216     break;
217     case 2:
218     val=-0.125*a1*c2;
219     break;
220     case 3:
221     val=-0.125*a1*b2;
222     break;
223     } break;
224     case 3:
225     switch (num_variable)
226     {
227     case 1:
228     val=0.125*b1*c2;
229     break;
230     case 2:
231     val=0.125*a1*c2;
232     break;
233     case 3:
234     val=-0.125*a1*b1;
235     break;
236     } break;
237     case 4:
238     switch (num_variable)
239     {
240     case 1:
241     val=-0.125*b1*c2;
242     break;
243     case 2:
244     val=0.125*a2*c2;
245     break;
246     case 3:
247     val=-0.125*a2*b1;
248     break;
249     } break;
250     case 5:
251     switch (num_variable)
252     {
253     case 1:
254     val=-0.125*b2*c1;
255     break;
256     case 2:
257     val=-0.125*a2*c1;
258     break;
259     case 3:
260     val=0.125*a2*b2;
261     break;
262     } break;
263     case 6:
264     switch (num_variable)
265     {
266     case 1:
267     val=0.125*b2*c1;
268     break;
269     case 2:
270     val=-0.125*a1*c1;
271     break;
272     case 3:
273     val=0.125*a1*b2;
274     break;
275     } break;
276     case 7:
277     switch (num_variable)
278     {
279     case 1:
280     val=0.125*b1*c1;
281     break;
282     case 2:
283     val=0.125*a1*c1;
284     break;
285     case 3:
286     val=0.125*a1*b1;
287     break;
288     } break;
289     case 8:
290     switch (num_variable)
291     {
292     case 1:
293     val=-0.125*b1*c1;
294     break;
295     case 2:
296     val=0.125*a2*c1;
297     break;
298     case 3:
299     val=0.125*a2*b1;
300     break;
301     } break;
302     }
303     return val;
304     }
305    
306     double FEM_HEXA8::get_jacobien(double* jac,double *uv,int& li,int& col,double unite)
307     {
308     li=3;col=3;
309     OT_TENSEUR Ni(3,8),Xi(8,3);
310     for (int i=0;i<3;i++)
311     for (int j=0;j<8;j++)
312     Ni(i,j)=get_fonction_derive_interpolation(j+1,i+1,uv);
313     for (int i=0;i<8;i++)
314     {
315     double *xyz=tab[i]->get_coord();
316     Xi(i,0)=xyz[0]*unite;
317     Xi(i,1)=xyz[1]*unite;
318     Xi(i,2)=xyz[2]*unite;
319     }
320     OT_TENSEUR jacobi=Ni*Xi;
321     jac[0]=jacobi(0,0).get_x();
322     jac[1]=jacobi(0,1).get_x();
323     jac[2]=jacobi(0,2).get_x();
324     jac[3]=jacobi(1,0).get_x();
325     jac[4]=jacobi(1,1).get_x();
326     jac[5]=jacobi(1,2).get_x();
327     jac[6]=jacobi(2,0).get_x();
328     jac[7]=jacobi(2,1).get_x();
329     jac[8]=jacobi(2,2).get_x();
330     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];
331     return det;
332    
333     }
334    
335     void FEM_HEXA8::get_inverse_jacob(double* j,double *uv,double unite)
336     {
337     double jac[9];
338     int li,col;
339     double detj=get_jacobien(jac,uv,li,col,unite);
340     j[0*3+0] = (jac[1*3+1]*jac[2*3+2]-jac[1*3+2]*jac[2*3+1])/(detj);
341     j[0*3+1] = -(jac[0*3+1]*jac[2*3+2]-jac[0*3+2]*jac[2*3+1])/(detj);
342     j[0*3+2] =-(-jac[0*3+1]*jac[1*3+2]+jac[0*3+2]*jac[1*3+1])/(detj);
343    
344     j[1*3+0] = -(jac[1*3+0]*jac[2*3+2]-jac[1*3+2]*jac[2*3+0])/(detj);
345     j[1*3+1] = (jac[0*3+0]*jac[2*3+2]-jac[0*3+2]*jac[2*3+0])/(detj);
346     j[1*3+2] = -(jac[0*3+0]*jac[1*3+2]-jac[0*3+2]*jac[1*3+0])/(detj);
347    
348     j[2*3+0] =-(-jac[1*3+0]*jac[2*3+1]+jac[1*3+1]*jac[2*3+0])/(detj);
349     j[2*3+1] = -(jac[0*3+0]*jac[2*3+1]-jac[0*3+1]*jac[2*3+0])/(detj);
350     j[2*3+2] = (jac[0*3+0]*jac[1*3+1]-jac[0*3+1]*jac[1*3+0])/(detj);
351    
352     }