ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_tetra4.cpp
Revision: 684
Committed: Tue Jun 30 15:21:22 2015 UTC (9 years, 10 months ago) by francois
File size: 15532 byte(s)
Log Message:
mise a jour de la boite englobante quand une coordonnees change  dans un fem_noeud

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