ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_tetra4.cpp
Revision: 638
Committed: Mon Jan 26 21:56:20 2015 UTC (10 years, 3 months ago) by francois
File size: 15455 byte(s)
Log Message:
ajout d'une méthode qui renvoit le degre max des fonctions d'interpolation et changement de nom de la fonction qui renvoie le nombre de fonction d'interpolation (ajout d'un get pour uniformiser)

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 francois 309 // fem_element34.cpp
16 francois 283 //
17     //------------------------------------------------------------
18     //------------------------------------------------------------
19     // COPYRIGHT 2000
20     // Version du 02/03/2006 � 11H22
21     //------------------------------------------------------------
22     //------------------------------------------------------------
23    
24    
25     #include "gestionversion.h"
26     #include "fem_tetra4.h"
27     #include "fem_noeud.h"
28     #include "mg_element_maillage.h"
29     #include "fem_maillage.h"
30 francois 405 #include "ot_mathematique.h"
31     #include "xfem_tetra4.h"
32     #include <math.h>
33 francois 283
34 francois 309 FEM_TETRA4::FEM_TETRA4(unsigned long num,class MG_ELEMENT_MAILLAGE* mai,class FEM_NOEUD** tabnoeud):FEM_ELEMENT3(num,mai),FEM_TEMPLATE_ELEMENT<4>(tabnoeud)
35 francois 283 {
36     if (liaison_topologique!=NULL)
37     if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
38 francois 309 tab[0]->get_lien_element3()->ajouter(this);
39     tab[1]->get_lien_element3()->ajouter(this);
40     tab[2]->get_lien_element3()->ajouter(this);
41     tab[3]->get_lien_element3()->ajouter(this);
42     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
43 francois 283 }
44    
45 francois 309 FEM_TETRA4::FEM_TETRA4(class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT3(mai),FEM_TEMPLATE_ELEMENT<4>(tabnoeud)
46 francois 283 {
47     if (liaison_topologique!=NULL)
48     if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
49 francois 309 tab[0]->get_lien_element3()->ajouter(this);
50     tab[1]->get_lien_element3()->ajouter(this);
51     tab[2]->get_lien_element3()->ajouter(this);
52     tab[3]->get_lien_element3()->ajouter(this);
53     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
54 francois 283 }
55 francois 378 FEM_TETRA4::FEM_TETRA4(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo,class FEM_NOEUD** tabnoeud):FEM_ELEMENT3(num,topo),FEM_TEMPLATE_ELEMENT<4>(tabnoeud)
56     {
57     if (liaison_topologique!=NULL)
58     if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
59     tab[0]->get_lien_element3()->ajouter(this);
60     tab[1]->get_lien_element3()->ajouter(this);
61     tab[2]->get_lien_element3()->ajouter(this);
62     tab[3]->get_lien_element3()->ajouter(this);
63     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
64     }
65 francois 283
66 francois 378 FEM_TETRA4::FEM_TETRA4(class MG_ELEMENT_TOPOLOGIQUE* topo,FEM_NOEUD** tabnoeud):FEM_ELEMENT3(topo),FEM_TEMPLATE_ELEMENT<4>(tabnoeud)
67     {
68     if (liaison_topologique!=NULL)
69     if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
70     tab[0]->get_lien_element3()->ajouter(this);
71     tab[1]->get_lien_element3()->ajouter(this);
72     tab[2]->get_lien_element3()->ajouter(this);
73     tab[3]->get_lien_element3()->ajouter(this);
74     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
75     }
76     FEM_TETRA4::FEM_TETRA4(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<4>(tabnoeud)
77     {
78     if (liaison_topologique!=NULL)
79     if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
80     tab[0]->get_lien_element3()->ajouter(this);
81     tab[1]->get_lien_element3()->ajouter(this);
82     tab[2]->get_lien_element3()->ajouter(this);
83     tab[3]->get_lien_element3()->ajouter(this);
84     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
85     }
86    
87     FEM_TETRA4::FEM_TETRA4(class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT3(topo,mai),FEM_TEMPLATE_ELEMENT<4>(tabnoeud)
88     {
89     if (liaison_topologique!=NULL)
90     if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
91     tab[0]->get_lien_element3()->ajouter(this);
92     tab[1]->get_lien_element3()->ajouter(this);
93     tab[2]->get_lien_element3()->ajouter(this);
94     tab[3]->get_lien_element3()->ajouter(this);
95     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
96     }
97    
98    
99    
100    
101    
102    
103    
104 francois 309 FEM_TETRA4::FEM_TETRA4(FEM_TETRA4& mdd):FEM_ELEMENT3(mdd),FEM_TEMPLATE_ELEMENT<4>(mdd)
105 francois 283 {
106     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
107 francois 309 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     get_fem_noeudpetitid()->get_lien_petit_element3()->ajouter(this);
112 francois 283 }
113    
114     FEM_TETRA4::~FEM_TETRA4()
115     {
116     if (liaison_topologique==NULL) return;
117     if (liaison_topologique->get_dimension()==0) liaison_topologique->get_lien_fem_maillage()->supprimer(this);
118 francois 309 tab[0]->get_lien_element3()->supprimer(this);
119     tab[1]->get_lien_element3()->supprimer(this);
120     tab[2]->get_lien_element3()->supprimer(this);
121     tab[3]->get_lien_element3()->supprimer(this);
122     get_fem_noeudpetitid()->get_lien_petit_element3()->supprimer(this);
123 francois 283 }
124    
125    
126     FEM_ELEMENT_MAILLAGE* FEM_TETRA4::dupliquer(FEM_MAILLAGE *femmai,long decalage)
127     {
128     FEM_NOEUD* tabnoeud[4];
129     tabnoeud[0]=femmai->get_fem_noeudid(tab[0]->get_id()+decalage);
130     tabnoeud[1]=femmai->get_fem_noeudid(tab[1]->get_id()+decalage);
131     tabnoeud[2]=femmai->get_fem_noeudid(tab[2]->get_id()+decalage);
132     tabnoeud[3]=femmai->get_fem_noeudid(tab[3]->get_id()+decalage);
133     FEM_TETRA4* tet=new FEM_TETRA4(get_id()+decalage,maillage,tabnoeud);
134 francois 309 femmai->ajouter_fem_element3(tet);
135 francois 283 return tet;
136     }
137    
138    
139     int FEM_TETRA4::get_type_entite(void)
140     {
141     return IDFEM_TETRA4;
142     }
143    
144     int FEM_TETRA4::get_dimension(void)
145     {
146     return 3;
147     }
148    
149    
150     int FEM_TETRA4::get_nb_fem_noeud(void)
151     {
152     return FEM_TEMPLATE_ELEMENT<4>::get_nb_fem_noeud();
153     }
154    
155     FEM_NOEUD* FEM_TETRA4::get_fem_noeud(int num)
156     {
157     return FEM_TEMPLATE_ELEMENT<4>::get_fem_noeud(num);
158     }
159    
160     void FEM_TETRA4::change_noeud(int num,FEM_NOEUD* noeud)
161     {
162     FEM_TEMPLATE_ELEMENT<4>::change_noeud(num,noeud);
163     }
164    
165    
166     BOITE_3D& FEM_TETRA4::get_boite_3D(void)
167     {
168     return FEM_TEMPLATE_ELEMENT<4>::get_boite_3D();
169     }
170 francois 637
171    
172     int FEM_TETRA4::get_nb_pt_gauss(int degre)
173     {
174     if (degre<2) return 1;
175     if (degre<3) return 4;
176     if (degre<4) return 5;
177     if (degre<6) return 15;
178     return 0;
179     }
180    
181     void FEM_TETRA4::get_pt_gauss(int degre,int num,double &w,double *uvw)
182     {
183     if (degre<2)
184     {
185     if (num==0) {w=0.166666666666667;uvw[0]=0.25;uvw[1]=0.25;uvw[2]=0.25;}
186     return;
187     }
188     if (degre<3)
189     {
190     if (num==0) {w=0.041666666666667;uvw[0]=0.138196601125011;uvw[1]=0.138196601125011;uvw[2]=0.138196601125011;}
191     if (num==1) {w=0.041666666666667;uvw[0]=0.138196601125011;uvw[1]=0.138196601125011;uvw[2]=0.585410196624968;}
192     if (num==2) {w=0.041666666666667;uvw[0]=0.138196601125011;uvw[1]=0.585410196624968;uvw[2]=0.138196601125011;}
193     if (num==3) {w=0.041666666666667;uvw[0]=0.585410196624968;uvw[1]=0.138196601125011;uvw[2]=0.138196601125011;}
194     return;
195     }
196     if (degre<4)
197     {
198     if (num==0) {w=-0.133333333333333;uvw[0]=0.25;uvw[1]=0.25;uvw[2]=0.25;}
199     if (num==1) {w=0.075;uvw[0]=0.166666666666667;uvw[1]=0.166666666666667;uvw[2]=0.166666666666667;}
200     if (num==2) {w=0.075;uvw[0]=0.166666666666667;uvw[1]=0.166666666666667;uvw[2]=0.5;}
201     if (num==3) {w=0.075;uvw[0]=0.166666666666667;uvw[1]=0.5;uvw[2]=0.166666666666667;}
202     if (num==4) {w=0.075;uvw[0]=0.5;uvw[1]=0.166666666666667;uvw[2]=0.166666666666667;}
203     return;
204     }
205     if (degre<6)
206     {
207     double a=0.25;
208     double b1=0.3197936278;
209     double b2=0.0919710781;
210     double c2=0.7240867658;
211     double c1=0.0406191165;
212     double d=0.0563508327;
213     double e=0.4436491673;
214     double w1=0.0197530864;
215     double w3=0.011989514;
216     double w2=0.0115113679;
217     double w4=0.0088183422;
218     if (num==0) {w=w1;uvw[0]=a;uvw[1]=a;uvw[2]=a;}
219     if (num==1) {w=w2;uvw[0]=b1;uvw[1]=b1;uvw[2]=b1;}
220     if (num==2) {w=w2;uvw[0]=c1;uvw[1]=b1;uvw[2]=b1;}
221     if (num==3) {w=w2;uvw[0]=b1;uvw[1]=c1;uvw[2]=b1;}
222     if (num==4) {w=w2;uvw[0]=b1;uvw[1]=b1;uvw[2]=c1;}
223     if (num==5) {w=w3;uvw[0]=b2;uvw[1]=b2;uvw[2]=b2;}
224     if (num==6) {w=w3;uvw[0]=c2;uvw[1]=b2;uvw[2]=b2;}
225     if (num==7) {w=w3;uvw[0]=b2;uvw[1]=c2;uvw[2]=b2;}
226     if (num==8) {w=w3;uvw[0]=b2;uvw[1]=b2;uvw[2]=c2;}
227     if (num==9) {w=w4;uvw[0]=e;uvw[1]=d;uvw[2]=d;}
228     if (num==10) {w=w4;uvw[0]=d;uvw[1]=e;uvw[2]=d;}
229     if (num==11) {w=w4;uvw[0]=d;uvw[1]=d;uvw[2]=e;}
230     if (num==12) {w=w4;uvw[0]=d;uvw[1]=e;uvw[2]=e;}
231     if (num==13) {w=w4;uvw[0]=e;uvw[1]=d;uvw[2]=e;}
232     if (num==14) {w=w4;uvw[0]=e;uvw[1]=e;uvw[2]=d;}
233     return;
234     }
235     }
236    
237    
238    
239 francois 638 int FEM_TETRA4::get_nb_fonction_interpolation(void)
240 francois 283 {
241     return 4;
242     }
243 francois 638 int FEM_TETRA4::get_degremax_fonction_interpolation(void)
244     {
245     return 1;
246     }
247 francois 283
248     double FEM_TETRA4::get_fonction_interpolation(int num,double *uv)
249     {
250     double val;
251     switch (num)
252     {
253     case 1:
254     val=1-uv[0]-uv[1]-uv[2];
255     break;
256     case 2:
257     val=uv[0];
258     break;
259     case 3:
260     val=uv[1];
261     break;
262     case 4:
263     val=uv[2];
264     break;
265     }
266     return val;
267     }
268    
269     double FEM_TETRA4::get_fonction_derive_interpolation(int num,int num_variable,double *uv)
270     {
271     double val;
272    
273     switch (num)
274     {
275     case 1:
276     switch (num_variable)
277     {
278     case 1:
279     val=-1;
280     break;
281     case 2:
282     val=-1;
283     break;
284     case 3:
285     val=-1;
286     break;
287     } break;
288     case 2:
289     switch (num_variable)
290     {
291     case 1:
292     val=1;
293     break;
294     case 2:
295     val=0;
296     break;
297     case 3:
298     val=0;
299     break;
300     }break;
301     case 3:
302     switch (num_variable)
303     {
304     case 1:
305     val=0;
306     break;
307     case 2:
308     val=1;
309     break;
310     case 3:
311     val=0;
312     break;
313     }break;
314     case 4:
315     switch (num_variable)
316     {
317     case 1:
318     val=0;
319     break;
320     case 2:
321     val=0;
322     break;
323     case 3:
324     val=1;
325     break;
326     }break;
327    
328     }
329     return val;
330     }
331    
332     double FEM_TETRA4::get_jacobien(double* jac,double *uv,int& li,int& col,double unite)
333     {
334     li=3;
335     col=3;
336     double *xyz1=tab[0]->get_coord();
337     double *xyz2=tab[1]->get_coord();
338     double *xyz3=tab[2]->get_coord();
339     double *xyz4=tab[3]->get_coord();
340    
341     jac[0]=unite*(xyz2[0] - xyz1[0]) ;
342     jac[1]=unite*(xyz2[1] - xyz1[1]) ;
343     jac[2]=unite*(xyz2[2] - xyz1[2]) ;
344    
345     jac[3]=unite*(xyz3[0] - xyz1[0]) ;
346     jac[4]=unite*(xyz3[1] - xyz1[1]) ;
347     jac[5]=unite*(xyz3[2] - xyz1[2]) ;
348    
349     jac[6]=unite*(xyz4[0] - xyz1[0]) ;
350     jac[7]=unite*(xyz4[1] - xyz1[1]) ;
351     jac[8]=unite*(xyz4[2] - xyz1[2]) ;
352    
353     double SIX_V= 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]);
354    
355     return SIX_V;
356     }
357    
358     void FEM_TETRA4::get_inverse_jacob(double* j,double *uv,double unite)
359     {
360     double *xyz1=tab[0]->get_coord();
361     double *xyz2=tab[1]->get_coord();
362     double *xyz3=tab[2]->get_coord();
363     double *xyz4=tab[3]->get_coord();
364    
365     double jac[9];
366    
367     jac[0*3+0]= unite*(xyz2[0] - xyz1[0]) ;
368     jac[0*3+1]= unite*(xyz2[1] - xyz1[1]) ;
369     jac[0*3+2]= unite*(xyz2[2] - xyz1[2]) ;
370    
371     jac[1*3+0]= unite*(xyz3[0] - xyz1[0]) ;
372     jac[1*3+1]= unite*(xyz3[1] - xyz1[1]) ;
373     jac[1*3+2]= unite*(xyz3[2] - xyz1[2]) ;
374    
375     jac[2*3+0]= unite*(xyz4[0] - xyz1[0]) ;
376     jac[2*3+1]= unite*(xyz4[1] - xyz1[1]) ;
377     jac[2*3+2]= unite*(xyz4[2] - xyz1[2]) ;
378    
379     double detj=(jac[0*3+0]*jac[1*3+1]*jac[2*3+2]-jac[0*3+0]*jac[1*3+2]*jac[2*3+1]-jac[1*3+0]*jac[0*3+1]*jac[2*3+2]+jac[1*3+0]*jac[0*3+2]*jac[2*3+1]+jac[2*3+0]*jac[0*3+1]*jac[1*3+2]-jac[2*3+0]*jac[0*3+2]*jac[1*3+1]);
380    
381     j[0*3+0] = (jac[1*3+1]*jac[2*3+2]-jac[1*3+2]*jac[2*3+1])/(detj);
382     j[0*3+1] = -(jac[0*3+1]*jac[2*3+2]-jac[0*3+2]*jac[2*3+1])/(detj);
383     j[0*3+2] =-(-jac[0*3+1]*jac[1*3+2]+jac[0*3+2]*jac[1*3+1])/(detj);
384    
385     j[1*3+0] = -(jac[1*3+0]*jac[2*3+2]-jac[1*3+2]*jac[2*3+0])/(detj);
386     j[1*3+1] = (jac[0*3+0]*jac[2*3+2]-jac[0*3+2]*jac[2*3+0])/(detj);
387     j[1*3+2] = -(jac[0*3+0]*jac[1*3+2]-jac[0*3+2]*jac[1*3+0])/(detj);
388    
389     j[2*3+0] =-(-jac[1*3+0]*jac[2*3+1]+jac[1*3+1]*jac[2*3+0])/(detj);
390     j[2*3+1] = -(jac[0*3+0]*jac[2*3+1]-jac[0*3+1]*jac[2*3+0])/(detj);
391     j[2*3+2] = (jac[0*3+0]*jac[1*3+1]-jac[0*3+1]*jac[1*3+0])/(detj);
392    
393    
394    
395     }
396    
397 francois 635 bool FEM_TETRA4::valide_parametre_EF(double* uvw)
398     {
399     if (uvw[0]>=-1e-10)
400     if (uvw[1]>=-1e-10)
401     if (uvw[2]>=-1e-10)
402     if (uvw[0]+uvw[1]+uvw[2]<=1.+1e-10)
403     return true;
404     return false;
405     }
406 francois 283 void FEM_TETRA4::enregistrer(std::ostream& o)
407     {
408 francois 378 if (maillage!=NULL)
409 francois 410 if (get_lien_topologie()!=NULL) o << "%" << get_id() << "=FEM_TETRA4($"<< get_lien_topologie()->get_id() << ",$" << maillage->get_id() << ",$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id() << ",$" << tab[3]->get_id() << ","<< get_etat(0) << ");" << std::endl;
410     else o << "%" << get_id() << "=FEM_TETRA4(NULL,$"<< maillage->get_id() << ",$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id() << ",$" << tab[3]->get_id() << ","<< get_etat(0) << ");" << std::endl;
411 francois 378 else
412 francois 410 if (get_lien_topologie()!=NULL) o << "%" << get_id() << "=FEM_TETRA4($"<< get_lien_topologie()->get_id() << ",NULL,$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id() << ",$" << tab[3]->get_id() << ","<< get_etat(0) << ");" << std::endl;
413     else o << "%" << get_id() << "=FEM_TETRA4(NULL,NULL,$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id() << ",$" << tab[3]->get_id() << ","<< get_etat(0) << ");" << std::endl;
414 francois 339
415 francois 283 }
416 francois 405
417    
418 francois 406 int FEM_TETRA4::verifie_validite_decoupage_xfem(double *vol)
419 francois 405 {
420     if (get_nb_xfem(3)==0) return 1;
421     FEM_NOEUD* no1=get_fem_noeud(0);
422     FEM_NOEUD* no2=get_fem_noeud(1);
423     FEM_NOEUD* no3=get_fem_noeud(2);
424     FEM_NOEUD* no4=get_fem_noeud(3);
425     OT_VECTEUR_3D vec1(no1->get_coord(),no2->get_coord());
426     OT_VECTEUR_3D vec2(no1->get_coord(),no3->get_coord());
427     OT_VECTEUR_3D vec3(no1->get_coord(),no4->get_coord());
428     double volume=(vec1&vec2)*vec3;
429 francois 407 int volumenul=0;
430 francois 405 for (int i=0;i<get_nb_xfem(3);i++)
431     {
432     XFEM_TETRA4* xtet=(XFEM_TETRA4*)get_xfem(3,i);
433     FEM_NOEUD* xno1=xtet->get_fem_noeud(0);
434     FEM_NOEUD* xno2=xtet->get_fem_noeud(1);
435     FEM_NOEUD* xno3=xtet->get_fem_noeud(2);
436     FEM_NOEUD* xno4=xtet->get_fem_noeud(3);
437     OT_VECTEUR_3D xvec1(xno1->get_coord(),xno2->get_coord());
438     OT_VECTEUR_3D xvec2(xno1->get_coord(),xno3->get_coord());
439     OT_VECTEUR_3D xvec3(xno1->get_coord(),xno4->get_coord());
440     double xvolume=(xvec1&xvec2)*xvec3;
441     volume=volume-xvolume;
442 francois 407 if (xvolume<1e-10*volume) volumenul=1;
443 francois 405 }
444 francois 406 if (vol!=NULL) *vol=volume;
445 francois 407 if (fabs(volume)>1e-12) return 2*volumenul;
446     return 1+2*volumenul;
447 francois 405 }