ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_tetra4.cpp
Revision: 407
Committed: Wed Jun 5 22:07:18 2013 UTC (11 years, 11 months ago) by francois
File size: 12691 byte(s)
Log Message: 
verification plus poussee des decoupage xfem

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     int FEM_TETRA4::nb_fonction_interpolation(void)
171     {
172     return 4;
173     }
174    
175    
176     double FEM_TETRA4::get_fonction_interpolation(int num,double *uv)
177     {
178     double val;
179     switch (num)
180     {
181     case 1:
182     val=1-uv[0]-uv[1]-uv[2];
183     break;
184     case 2:
185     val=uv[0];
186     break;
187     case 3:
188     val=uv[1];
189     break;
190     case 4:
191     val=uv[2];
192     break;
193     }
194     return val;
195     }
196    
197     double FEM_TETRA4::get_fonction_derive_interpolation(int num,int num_variable,double *uv)
198     {
199     double val;
200    
201     switch (num)
202     {
203     case 1:
204     switch (num_variable)
205     {
206     case 1:
207     val=-1;
208     break;
209     case 2:
210     val=-1;
211     break;
212     case 3:
213     val=-1;
214     break;
215     } break;
216     case 2:
217     switch (num_variable)
218     {
219     case 1:
220     val=1;
221     break;
222     case 2:
223     val=0;
224     break;
225     case 3:
226     val=0;
227     break;
228     }break;
229     case 3:
230     switch (num_variable)
231     {
232     case 1:
233     val=0;
234     break;
235     case 2:
236     val=1;
237     break;
238     case 3:
239     val=0;
240     break;
241     }break;
242     case 4:
243     switch (num_variable)
244     {
245     case 1:
246     val=0;
247     break;
248     case 2:
249     val=0;
250     break;
251     case 3:
252     val=1;
253     break;
254     }break;
255    
256     }
257     return val;
258     }
259    
260     double FEM_TETRA4::get_jacobien(double* jac,double *uv,int& li,int& col,double unite)
261     {
262     li=3;
263     col=3;
264     double *xyz1=tab[0]->get_coord();
265     double *xyz2=tab[1]->get_coord();
266     double *xyz3=tab[2]->get_coord();
267     double *xyz4=tab[3]->get_coord();
268    
269     jac[0]=unite*(xyz2[0] - xyz1[0]) ;
270     jac[1]=unite*(xyz2[1] - xyz1[1]) ;
271     jac[2]=unite*(xyz2[2] - xyz1[2]) ;
272    
273     jac[3]=unite*(xyz3[0] - xyz1[0]) ;
274     jac[4]=unite*(xyz3[1] - xyz1[1]) ;
275     jac[5]=unite*(xyz3[2] - xyz1[2]) ;
276    
277     jac[6]=unite*(xyz4[0] - xyz1[0]) ;
278     jac[7]=unite*(xyz4[1] - xyz1[1]) ;
279     jac[8]=unite*(xyz4[2] - xyz1[2]) ;
280    
281     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]);
282    
283     return SIX_V;
284     }
285    
286     void FEM_TETRA4::get_inverse_jacob(double* j,double *uv,double unite)
287     {
288     double *xyz1=tab[0]->get_coord();
289     double *xyz2=tab[1]->get_coord();
290     double *xyz3=tab[2]->get_coord();
291     double *xyz4=tab[3]->get_coord();
292    
293     double jac[9];
294    
295     jac[0*3+0]= unite*(xyz2[0] - xyz1[0]) ;
296     jac[0*3+1]= unite*(xyz2[1] - xyz1[1]) ;
297     jac[0*3+2]= unite*(xyz2[2] - xyz1[2]) ;
298    
299     jac[1*3+0]= unite*(xyz3[0] - xyz1[0]) ;
300     jac[1*3+1]= unite*(xyz3[1] - xyz1[1]) ;
301     jac[1*3+2]= unite*(xyz3[2] - xyz1[2]) ;
302    
303     jac[2*3+0]= unite*(xyz4[0] - xyz1[0]) ;
304     jac[2*3+1]= unite*(xyz4[1] - xyz1[1]) ;
305     jac[2*3+2]= unite*(xyz4[2] - xyz1[2]) ;
306    
307     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]);
308    
309     j[0*3+0] = (jac[1*3+1]*jac[2*3+2]-jac[1*3+2]*jac[2*3+1])/(detj);
310     j[0*3+1] = -(jac[0*3+1]*jac[2*3+2]-jac[0*3+2]*jac[2*3+1])/(detj);
311     j[0*3+2] =-(-jac[0*3+1]*jac[1*3+2]+jac[0*3+2]*jac[1*3+1])/(detj);
312    
313     j[1*3+0] = -(jac[1*3+0]*jac[2*3+2]-jac[1*3+2]*jac[2*3+0])/(detj);
314     j[1*3+1] = (jac[0*3+0]*jac[2*3+2]-jac[0*3+2]*jac[2*3+0])/(detj);
315     j[1*3+2] = -(jac[0*3+0]*jac[1*3+2]-jac[0*3+2]*jac[1*3+0])/(detj);
316    
317     j[2*3+0] =-(-jac[1*3+0]*jac[2*3+1]+jac[1*3+1]*jac[2*3+0])/(detj);
318     j[2*3+1] = -(jac[0*3+0]*jac[2*3+1]-jac[0*3+1]*jac[2*3+0])/(detj);
319     j[2*3+2] = (jac[0*3+0]*jac[1*3+1]-jac[0*3+1]*jac[1*3+0])/(detj);
320    
321    
322    
323     }
324    
325    
326     void FEM_TETRA4::enregistrer(std::ostream& o)
327     {
328 francois 378 if (maillage!=NULL)
329     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() << ");" << std::endl;
330     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() << ");" << std::endl;
331     else
332     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() << ");" << std::endl;
333     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() << ");" << std::endl;
334 francois 339
335 francois 283 }
336 francois 405
337    
338 francois 406 int FEM_TETRA4::verifie_validite_decoupage_xfem(double *vol)
339 francois 405 {
340     if (get_nb_xfem(3)==0) return 1;
341     FEM_NOEUD* no1=get_fem_noeud(0);
342     FEM_NOEUD* no2=get_fem_noeud(1);
343     FEM_NOEUD* no3=get_fem_noeud(2);
344     FEM_NOEUD* no4=get_fem_noeud(3);
345     OT_VECTEUR_3D vec1(no1->get_coord(),no2->get_coord());
346     OT_VECTEUR_3D vec2(no1->get_coord(),no3->get_coord());
347     OT_VECTEUR_3D vec3(no1->get_coord(),no4->get_coord());
348     double volume=(vec1&vec2)*vec3;
349 francois 407 int volumenul=0;
350 francois 405 for (int i=0;i<get_nb_xfem(3);i++)
351     {
352     XFEM_TETRA4* xtet=(XFEM_TETRA4*)get_xfem(3,i);
353     FEM_NOEUD* xno1=xtet->get_fem_noeud(0);
354     FEM_NOEUD* xno2=xtet->get_fem_noeud(1);
355     FEM_NOEUD* xno3=xtet->get_fem_noeud(2);
356     FEM_NOEUD* xno4=xtet->get_fem_noeud(3);
357     OT_VECTEUR_3D xvec1(xno1->get_coord(),xno2->get_coord());
358     OT_VECTEUR_3D xvec2(xno1->get_coord(),xno3->get_coord());
359     OT_VECTEUR_3D xvec3(xno1->get_coord(),xno4->get_coord());
360     double xvolume=(xvec1&xvec2)*xvec3;
361     volume=volume-xvolume;
362 francois 407 if (xvolume<1e-10*volume) volumenul=1;
363 francois 405 }
364 francois 406 if (vol!=NULL) *vol=volume;
365 francois 407 if (fabs(volume)>1e-12) return 2*volumenul;
366     return 1+2*volumenul;
367 francois 405 }