ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_hexa8.cpp
Revision: 754
Committed: Thu Oct 29 21:49:15 2015 UTC (9 years, 6 months ago) by francois
File size: 15326 byte(s)
Log Message:
Centralisation et unification de la définition des points de gauss des elements finis et arrangement dans le meme ordre que 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    
166    
167 francois 635 bool FEM_HEXA8::valide_parametre_EF(double* uvw)
168     {
169 francois 674 if (uvw[0]>=-1.-1e-10)
170     if (uvw[1]>=-1.-1e-10)
171     if (uvw[2]>=-1.-1e-10)
172 francois 635 if (uvw[0]<=1.+1e-10)
173     if (uvw[1]<=1.+1e-10)
174     if (uvw[2]<=1.+1e-10)
175     return true;
176     return false;
177     }
178 francois 310
179 francois 684 void FEM_HEXA8::reinit_boite_3D(void)
180     {
181     FEM_TEMPLATE_ELEMENT<8>::reinit_boite_3D();
182     }
183 francois 310
184     int FEM_HEXA8::get_type_entite(void)
185     {
186     return IDFEM_HEXA8;
187     }
188    
189     int FEM_HEXA8::get_dimension(void)
190     {
191     return 3;
192     }
193    
194    
195    
196     int FEM_HEXA8::get_nb_fem_noeud(void)
197     {
198     return FEM_TEMPLATE_ELEMENT<8>::get_nb_fem_noeud();
199     }
200    
201     FEM_NOEUD* FEM_HEXA8::get_fem_noeud(int num)
202     {
203     return FEM_TEMPLATE_ELEMENT<8>::get_fem_noeud(num);
204     }
205    
206     void FEM_HEXA8::change_noeud(int num,FEM_NOEUD* noeud)
207     {
208     FEM_TEMPLATE_ELEMENT<8>::change_noeud(num,noeud);
209     }
210    
211     BOITE_3D& FEM_HEXA8::get_boite_3D(void)
212     {
213     return FEM_TEMPLATE_ELEMENT<8>::get_boite_3D();
214     }
215    
216    
217     void FEM_HEXA8::enregistrer(std::ostream& o)
218     {
219 francois 378 if (maillage!=NULL)
220     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;
221     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;
222     else
223     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;
224     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;
225    
226 francois 310 }
227 francois 637
228    
229     int FEM_HEXA8::get_nb_pt_gauss(int degre)
230     {
231 francois 754 return OT_POINTS_GAUSS::get_nb_point_hexa_prod(degre);
232 francois 637 }
233    
234     void FEM_HEXA8::get_pt_gauss(int degre,int num,double &w,double *uvw)
235     {
236 francois 754 return OT_POINTS_GAUSS::get_pt_gauss_hex_prod(degre,num,w,uvw);
237 francois 637 }
238 francois 754
239 francois 638 int FEM_HEXA8::get_nb_fonction_interpolation(void)
240 francois 310 {
241     return 8;
242     }
243    
244 francois 638 int FEM_HEXA8::get_degremax_fonction_interpolation(void)
245     {
246     return 3;
247     }
248 francois 310 double FEM_HEXA8::get_fonction_interpolation(int num,double *uv)
249     {
250     double val;
251     double a1=1.+uv[0],a2=1.-uv[0];
252     double b1=1.+uv[1],b2=1.-uv[1];
253     double c1=1.+uv[2],c2=1.-uv[2];
254     switch (num)
255     {
256     case 1:
257     val=0.125*a2*b2*c2;
258     break;
259     case 2:
260     val=0.125*a1*b2*c2;
261     break;
262     case 3:
263     val=0.125*a1*b1*c2;
264     break;
265     case 4:
266     val=0.125*a2*b1*c2;
267     break;
268     case 5:
269     val=0.125*a2*b2*c1;
270     break;
271     case 6:
272     val=0.125*a1*b2*c1;
273     break;
274     case 7:
275     val=0.125*a1*b1*c1;
276     break;
277     case 8:
278     val=0.125*a2*b1*c1;
279     break;
280     }
281     return val;
282     }
283    
284     double FEM_HEXA8::get_fonction_derive_interpolation(int num,int num_variable,double *uv)
285     {
286     double val;
287     double a1=1.+uv[0],a2=1.-uv[0];
288     double b1=1.+uv[1],b2=1.-uv[1];
289     double c1=1.+uv[2],c2=1.-uv[2];
290     switch (num)
291     {
292     case 1:
293     switch (num_variable)
294     {
295     case 1:
296     val=-0.125*b2*c2;
297     break;
298     case 2:
299     val=-0.125*a2*c2;
300     break;
301     case 3:
302     val=-0.125*a2*b2;
303     break;
304     } break;
305     case 2:
306     switch (num_variable)
307     {
308     case 1:
309     val=0.125*b2*c2;
310     break;
311     case 2:
312     val=-0.125*a1*c2;
313     break;
314     case 3:
315     val=-0.125*a1*b2;
316     break;
317     } break;
318     case 3:
319     switch (num_variable)
320     {
321     case 1:
322     val=0.125*b1*c2;
323     break;
324     case 2:
325     val=0.125*a1*c2;
326     break;
327     case 3:
328     val=-0.125*a1*b1;
329     break;
330     } break;
331     case 4:
332     switch (num_variable)
333     {
334     case 1:
335     val=-0.125*b1*c2;
336     break;
337     case 2:
338     val=0.125*a2*c2;
339     break;
340     case 3:
341     val=-0.125*a2*b1;
342     break;
343     } break;
344     case 5:
345     switch (num_variable)
346     {
347     case 1:
348     val=-0.125*b2*c1;
349     break;
350     case 2:
351     val=-0.125*a2*c1;
352     break;
353     case 3:
354     val=0.125*a2*b2;
355     break;
356     } break;
357     case 6:
358     switch (num_variable)
359     {
360     case 1:
361     val=0.125*b2*c1;
362     break;
363     case 2:
364     val=-0.125*a1*c1;
365     break;
366     case 3:
367     val=0.125*a1*b2;
368     break;
369     } break;
370     case 7:
371     switch (num_variable)
372     {
373     case 1:
374     val=0.125*b1*c1;
375     break;
376     case 2:
377     val=0.125*a1*c1;
378     break;
379     case 3:
380     val=0.125*a1*b1;
381     break;
382     } break;
383     case 8:
384     switch (num_variable)
385     {
386     case 1:
387     val=-0.125*b1*c1;
388     break;
389     case 2:
390     val=0.125*a2*c1;
391     break;
392     case 3:
393     val=0.125*a2*b1;
394     break;
395     } break;
396     }
397     return val;
398     }
399    
400     double FEM_HEXA8::get_jacobien(double* jac,double *uv,int& li,int& col,double unite)
401     {
402     li=3;col=3;
403     OT_TENSEUR Ni(3,8),Xi(8,3);
404     for (int i=0;i<3;i++)
405     for (int j=0;j<8;j++)
406     Ni(i,j)=get_fonction_derive_interpolation(j+1,i+1,uv);
407     for (int i=0;i<8;i++)
408     {
409     double *xyz=tab[i]->get_coord();
410     Xi(i,0)=xyz[0]*unite;
411     Xi(i,1)=xyz[1]*unite;
412     Xi(i,2)=xyz[2]*unite;
413     }
414     OT_TENSEUR jacobi=Ni*Xi;
415     jac[0]=jacobi(0,0).get_x();
416     jac[1]=jacobi(0,1).get_x();
417     jac[2]=jacobi(0,2).get_x();
418     jac[3]=jacobi(1,0).get_x();
419     jac[4]=jacobi(1,1).get_x();
420     jac[5]=jacobi(1,2).get_x();
421     jac[6]=jacobi(2,0).get_x();
422     jac[7]=jacobi(2,1).get_x();
423     jac[8]=jacobi(2,2).get_x();
424     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];
425     return det;
426    
427     }
428    
429     void FEM_HEXA8::get_inverse_jacob(double* j,double *uv,double unite)
430     {
431     double jac[9];
432     int li,col;
433     double detj=get_jacobien(jac,uv,li,col,unite);
434     j[0*3+0] = (jac[1*3+1]*jac[2*3+2]-jac[1*3+2]*jac[2*3+1])/(detj);
435     j[0*3+1] = -(jac[0*3+1]*jac[2*3+2]-jac[0*3+2]*jac[2*3+1])/(detj);
436     j[0*3+2] =-(-jac[0*3+1]*jac[1*3+2]+jac[0*3+2]*jac[1*3+1])/(detj);
437    
438     j[1*3+0] = -(jac[1*3+0]*jac[2*3+2]-jac[1*3+2]*jac[2*3+0])/(detj);
439     j[1*3+1] = (jac[0*3+0]*jac[2*3+2]-jac[0*3+2]*jac[2*3+0])/(detj);
440     j[1*3+2] = -(jac[0*3+0]*jac[1*3+2]-jac[0*3+2]*jac[1*3+0])/(detj);
441    
442     j[2*3+0] =-(-jac[1*3+0]*jac[2*3+1]+jac[1*3+1]*jac[2*3+0])/(detj);
443     j[2*3+1] = -(jac[0*3+0]*jac[2*3+1]-jac[0*3+1]*jac[2*3+0])/(detj);
444     j[2*3+2] = (jac[0*3+0]*jac[1*3+1]-jac[0*3+1]*jac[1*3+0])/(detj);
445    
446 francois 405 }
447    
448    
449    
450 francois 406 int FEM_HEXA8::verifie_validite_decoupage_xfem(double *vol)
451 francois 405 {
452    
453 francois 310 }