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, 7 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

# 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::get_nb_fonction_interpolation(void)
240 {
241 return 4;
242 }
243 int FEM_TETRA4::get_degremax_fonction_interpolation(void)
244 {
245 return 1;
246 }
247
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 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 void FEM_TETRA4::enregistrer(std::ostream& o)
407 {
408 if (maillage!=NULL)
409 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 else
412 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
415 }
416
417
418 int FEM_TETRA4::verifie_validite_decoupage_xfem(double *vol)
419 {
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 int volumenul=0;
430 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 if (xvolume<1e-10*volume) volumenul=1;
443 }
444 if (vol!=NULL) *vol=volume;
445 if (fabs(volume)>1e-12) return 2*volumenul;
446 return 1+2*volumenul;
447 }