ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_tetra4.cpp
Revision: 406
Committed: Wed Jun 5 20:13:26 2013 UTC (12 years, 2 months ago) by francois
File size: 12609 byte(s)
Log Message:
Procedure d'autoevaluation du decoupage xfem

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 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 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
335 }
336
337
338 int FEM_TETRA4::verifie_validite_decoupage_xfem(double *vol)
339 {
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 for (int i=0;i<get_nb_xfem(3);i++)
350 {
351 XFEM_TETRA4* xtet=(XFEM_TETRA4*)get_xfem(3,i);
352 FEM_NOEUD* xno1=xtet->get_fem_noeud(0);
353 FEM_NOEUD* xno2=xtet->get_fem_noeud(1);
354 FEM_NOEUD* xno3=xtet->get_fem_noeud(2);
355 FEM_NOEUD* xno4=xtet->get_fem_noeud(3);
356 OT_VECTEUR_3D xvec1(xno1->get_coord(),xno2->get_coord());
357 OT_VECTEUR_3D xvec2(xno1->get_coord(),xno3->get_coord());
358 OT_VECTEUR_3D xvec3(xno1->get_coord(),xno4->get_coord());
359 double xvolume=(xvec1&xvec2)*xvec3;
360 volume=volume-xvolume;
361 }
362 if (vol!=NULL) *vol=volume;
363 if (fabs(volume)>1e-12) return 0;
364 return 1;
365 }