ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_tetra4.cpp
Revision: 637
Committed: Mon Jan 26 21:26:37 2015 UTC (10 years, 7 months ago) by francois
File size: 15373 byte(s)
Log Message:
entree des points de gauss pour les différents éléments finis depuis reference JCC

File Contents

# Content
1 //------------------------------------------------------------
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_element34.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_tetra4.h"
27 #include "fem_noeud.h"
28 #include "mg_element_maillage.h"
29 #include "fem_maillage.h"
30 #include "ot_mathematique.h"
31 #include "xfem_tetra4.h"
32 #include <math.h>
33
34 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 {
36 if (liaison_topologique!=NULL)
37 if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
38 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 }
44
45 FEM_TETRA4::FEM_TETRA4(class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT3(mai),FEM_TEMPLATE_ELEMENT<4>(tabnoeud)
46 {
47 if (liaison_topologique!=NULL)
48 if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
49 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 }
55 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
66 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 FEM_TETRA4::FEM_TETRA4(FEM_TETRA4& mdd):FEM_ELEMENT3(mdd),FEM_TEMPLATE_ELEMENT<4>(mdd)
105 {
106 if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==3) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
107 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 }
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 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 }
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 femmai->ajouter_fem_element3(tet);
135 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
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 int FEM_TETRA4::nb_fonction_interpolation(void)
240 {
241 return 4;
242 }
243
244
245 double FEM_TETRA4::get_fonction_interpolation(int num,double *uv)
246 {
247 double val;
248 switch (num)
249 {
250 case 1:
251 val=1-uv[0]-uv[1]-uv[2];
252 break;
253 case 2:
254 val=uv[0];
255 break;
256 case 3:
257 val=uv[1];
258 break;
259 case 4:
260 val=uv[2];
261 break;
262 }
263 return val;
264 }
265
266 double FEM_TETRA4::get_fonction_derive_interpolation(int num,int num_variable,double *uv)
267 {
268 double val;
269
270 switch (num)
271 {
272 case 1:
273 switch (num_variable)
274 {
275 case 1:
276 val=-1;
277 break;
278 case 2:
279 val=-1;
280 break;
281 case 3:
282 val=-1;
283 break;
284 } break;
285 case 2:
286 switch (num_variable)
287 {
288 case 1:
289 val=1;
290 break;
291 case 2:
292 val=0;
293 break;
294 case 3:
295 val=0;
296 break;
297 }break;
298 case 3:
299 switch (num_variable)
300 {
301 case 1:
302 val=0;
303 break;
304 case 2:
305 val=1;
306 break;
307 case 3:
308 val=0;
309 break;
310 }break;
311 case 4:
312 switch (num_variable)
313 {
314 case 1:
315 val=0;
316 break;
317 case 2:
318 val=0;
319 break;
320 case 3:
321 val=1;
322 break;
323 }break;
324
325 }
326 return val;
327 }
328
329 double FEM_TETRA4::get_jacobien(double* jac,double *uv,int& li,int& col,double unite)
330 {
331 li=3;
332 col=3;
333 double *xyz1=tab[0]->get_coord();
334 double *xyz2=tab[1]->get_coord();
335 double *xyz3=tab[2]->get_coord();
336 double *xyz4=tab[3]->get_coord();
337
338 jac[0]=unite*(xyz2[0] - xyz1[0]) ;
339 jac[1]=unite*(xyz2[1] - xyz1[1]) ;
340 jac[2]=unite*(xyz2[2] - xyz1[2]) ;
341
342 jac[3]=unite*(xyz3[0] - xyz1[0]) ;
343 jac[4]=unite*(xyz3[1] - xyz1[1]) ;
344 jac[5]=unite*(xyz3[2] - xyz1[2]) ;
345
346 jac[6]=unite*(xyz4[0] - xyz1[0]) ;
347 jac[7]=unite*(xyz4[1] - xyz1[1]) ;
348 jac[8]=unite*(xyz4[2] - xyz1[2]) ;
349
350 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]);
351
352 return SIX_V;
353 }
354
355 void FEM_TETRA4::get_inverse_jacob(double* j,double *uv,double unite)
356 {
357 double *xyz1=tab[0]->get_coord();
358 double *xyz2=tab[1]->get_coord();
359 double *xyz3=tab[2]->get_coord();
360 double *xyz4=tab[3]->get_coord();
361
362 double jac[9];
363
364 jac[0*3+0]= unite*(xyz2[0] - xyz1[0]) ;
365 jac[0*3+1]= unite*(xyz2[1] - xyz1[1]) ;
366 jac[0*3+2]= unite*(xyz2[2] - xyz1[2]) ;
367
368 jac[1*3+0]= unite*(xyz3[0] - xyz1[0]) ;
369 jac[1*3+1]= unite*(xyz3[1] - xyz1[1]) ;
370 jac[1*3+2]= unite*(xyz3[2] - xyz1[2]) ;
371
372 jac[2*3+0]= unite*(xyz4[0] - xyz1[0]) ;
373 jac[2*3+1]= unite*(xyz4[1] - xyz1[1]) ;
374 jac[2*3+2]= unite*(xyz4[2] - xyz1[2]) ;
375
376 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]);
377
378 j[0*3+0] = (jac[1*3+1]*jac[2*3+2]-jac[1*3+2]*jac[2*3+1])/(detj);
379 j[0*3+1] = -(jac[0*3+1]*jac[2*3+2]-jac[0*3+2]*jac[2*3+1])/(detj);
380 j[0*3+2] =-(-jac[0*3+1]*jac[1*3+2]+jac[0*3+2]*jac[1*3+1])/(detj);
381
382 j[1*3+0] = -(jac[1*3+0]*jac[2*3+2]-jac[1*3+2]*jac[2*3+0])/(detj);
383 j[1*3+1] = (jac[0*3+0]*jac[2*3+2]-jac[0*3+2]*jac[2*3+0])/(detj);
384 j[1*3+2] = -(jac[0*3+0]*jac[1*3+2]-jac[0*3+2]*jac[1*3+0])/(detj);
385
386 j[2*3+0] =-(-jac[1*3+0]*jac[2*3+1]+jac[1*3+1]*jac[2*3+0])/(detj);
387 j[2*3+1] = -(jac[0*3+0]*jac[2*3+1]-jac[0*3+1]*jac[2*3+0])/(detj);
388 j[2*3+2] = (jac[0*3+0]*jac[1*3+1]-jac[0*3+1]*jac[1*3+0])/(detj);
389
390
391
392 }
393
394 bool FEM_TETRA4::valide_parametre_EF(double* uvw)
395 {
396 if (uvw[0]>=-1e-10)
397 if (uvw[1]>=-1e-10)
398 if (uvw[2]>=-1e-10)
399 if (uvw[0]+uvw[1]+uvw[2]<=1.+1e-10)
400 return true;
401 return false;
402 }
403 void FEM_TETRA4::enregistrer(std::ostream& o)
404 {
405 if (maillage!=NULL)
406 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;
407 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;
408 else
409 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;
410 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;
411
412 }
413
414
415 int FEM_TETRA4::verifie_validite_decoupage_xfem(double *vol)
416 {
417 if (get_nb_xfem(3)==0) return 1;
418 FEM_NOEUD* no1=get_fem_noeud(0);
419 FEM_NOEUD* no2=get_fem_noeud(1);
420 FEM_NOEUD* no3=get_fem_noeud(2);
421 FEM_NOEUD* no4=get_fem_noeud(3);
422 OT_VECTEUR_3D vec1(no1->get_coord(),no2->get_coord());
423 OT_VECTEUR_3D vec2(no1->get_coord(),no3->get_coord());
424 OT_VECTEUR_3D vec3(no1->get_coord(),no4->get_coord());
425 double volume=(vec1&vec2)*vec3;
426 int volumenul=0;
427 for (int i=0;i<get_nb_xfem(3);i++)
428 {
429 XFEM_TETRA4* xtet=(XFEM_TETRA4*)get_xfem(3,i);
430 FEM_NOEUD* xno1=xtet->get_fem_noeud(0);
431 FEM_NOEUD* xno2=xtet->get_fem_noeud(1);
432 FEM_NOEUD* xno3=xtet->get_fem_noeud(2);
433 FEM_NOEUD* xno4=xtet->get_fem_noeud(3);
434 OT_VECTEUR_3D xvec1(xno1->get_coord(),xno2->get_coord());
435 OT_VECTEUR_3D xvec2(xno1->get_coord(),xno3->get_coord());
436 OT_VECTEUR_3D xvec3(xno1->get_coord(),xno4->get_coord());
437 double xvolume=(xvec1&xvec2)*xvec3;
438 volume=volume-xvolume;
439 if (xvolume<1e-10*volume) volumenul=1;
440 }
441 if (vol!=NULL) *vol=volume;
442 if (fabs(volume)>1e-12) return 2*volumenul;
443 return 1+2*volumenul;
444 }