ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_triangle3.cpp
Revision: 647
Committed: Fri Feb 6 23:18:47 2015 UTC (10 years, 3 months ago) by francois
File size: 12371 byte(s)
Log Message: 
Amelioration du calcul du jacobien et inverse dans FEM_TRIANGLE3

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     // fem_triangle3.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_triangle3.h"
27     #include "fem_maillage.h"
28     #include "fem_noeud.h"
29     #include "mg_element_maillage.h"
30 gervaislavoie 382 #include "math.h"
31 francois 283
32    
33 francois 309 FEM_TRIANGLE3::FEM_TRIANGLE3(unsigned long num,class MG_ELEMENT_MAILLAGE* mai,class FEM_NOEUD** tabnoeud):FEM_ELEMENT2(num,mai),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
34 francois 283 {
35     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
36 francois 309 tab[0]->get_lien_element2()->ajouter(this);
37     tab[1]->get_lien_element2()->ajouter(this);
38     tab[2]->get_lien_element2()->ajouter(this);
39     get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
40 francois 283 }
41    
42 francois 309 FEM_TRIANGLE3::FEM_TRIANGLE3(class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(mai),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
43 francois 283 {
44     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
45 francois 309 tab[0]->get_lien_element2()->ajouter(this);
46     tab[1]->get_lien_element2()->ajouter(this);
47     tab[2]->get_lien_element2()->ajouter(this);
48     get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
49 francois 283 }
50 francois 378 FEM_TRIANGLE3::FEM_TRIANGLE3(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo,class FEM_NOEUD** tabnoeud):FEM_ELEMENT2(num,topo),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
51     {
52     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
53     tab[0]->get_lien_element2()->ajouter(this);
54     tab[1]->get_lien_element2()->ajouter(this);
55     tab[2]->get_lien_element2()->ajouter(this);
56     get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
57     }
58 francois 283
59 francois 378 FEM_TRIANGLE3::FEM_TRIANGLE3(class MG_ELEMENT_TOPOLOGIQUE* topo,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(topo),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
60     {
61     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
62     tab[0]->get_lien_element2()->ajouter(this);
63     tab[1]->get_lien_element2()->ajouter(this);
64     tab[2]->get_lien_element2()->ajouter(this);
65     get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
66     }
67     FEM_TRIANGLE3::FEM_TRIANGLE3(unsigned long num,class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai,class FEM_NOEUD** tabnoeud):FEM_ELEMENT2(num,topo,mai),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
68     {
69     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
70     tab[0]->get_lien_element2()->ajouter(this);
71     tab[1]->get_lien_element2()->ajouter(this);
72     tab[2]->get_lien_element2()->ajouter(this);
73     get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
74     }
75    
76     FEM_TRIANGLE3::FEM_TRIANGLE3(class MG_ELEMENT_TOPOLOGIQUE* topo,class MG_ELEMENT_MAILLAGE* mai,FEM_NOEUD** tabnoeud):FEM_ELEMENT2(topo,mai),FEM_TEMPLATE_ELEMENT<3>(tabnoeud)
77     {
78     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
79     tab[0]->get_lien_element2()->ajouter(this);
80     tab[1]->get_lien_element2()->ajouter(this);
81     tab[2]->get_lien_element2()->ajouter(this);
82     get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
83     }
84 francois 309 FEM_TRIANGLE3::FEM_TRIANGLE3(FEM_TRIANGLE3& mdd):FEM_ELEMENT2(mdd),FEM_TEMPLATE_ELEMENT<3>(mdd)
85 francois 283 {
86     if (liaison_topologique!=NULL) if (liaison_topologique->get_dimension()==2) liaison_topologique->get_lien_fem_maillage()->ajouter(this);
87 francois 309 tab[0]->get_lien_element2()->ajouter(this);
88     tab[1]->get_lien_element2()->ajouter(this);
89     tab[2]->get_lien_element2()->ajouter(this);
90     get_fem_noeudpetitid()->get_lien_petit_element2()->ajouter(this);
91 francois 283 }
92     FEM_TRIANGLE3::~FEM_TRIANGLE3()
93     {
94     if (liaison_topologique==NULL) return;
95     if (liaison_topologique->get_dimension()==0) liaison_topologique->get_lien_fem_maillage()->supprimer(this);
96 francois 309 tab[0]->get_lien_element2()->supprimer(this);
97     tab[1]->get_lien_element2()->supprimer(this);
98     tab[2]->get_lien_element2()->supprimer(this);
99     get_fem_noeudpetitid()->get_lien_petit_element2()->supprimer(this);
100 francois 283 }
101    
102    
103     FEM_ELEMENT_MAILLAGE* FEM_TRIANGLE3::dupliquer(FEM_MAILLAGE *femmai,long decalage)
104     {
105     FEM_NOEUD* tabnoeud[3];
106     tabnoeud[0]=femmai->get_fem_noeudid(tab[0]->get_id()+decalage);
107     tabnoeud[1]=femmai->get_fem_noeudid(tab[1]->get_id()+decalage);
108     tabnoeud[2]=femmai->get_fem_noeudid(tab[2]->get_id()+decalage);
109     FEM_TRIANGLE3* tri=new FEM_TRIANGLE3(get_id()+decalage,maillage,tabnoeud);
110 francois 309 femmai->ajouter_fem_element2(tri);
111 francois 283 return tri;
112     }
113    
114    
115    
116    
117    
118     int FEM_TRIANGLE3::get_type_entite(void)
119     {
120     return IDFEM_TRIANGLE3;
121     }
122    
123     int FEM_TRIANGLE3::get_dimension(void)
124     {
125     return 2;
126     }
127    
128    
129    
130     int FEM_TRIANGLE3::get_nb_fem_noeud(void)
131     {
132     return FEM_TEMPLATE_ELEMENT<3>::get_nb_fem_noeud();
133     }
134    
135     FEM_NOEUD* FEM_TRIANGLE3::get_fem_noeud(int num)
136     {
137     return FEM_TEMPLATE_ELEMENT<3>::get_fem_noeud(num);
138     }
139    
140     void FEM_TRIANGLE3::change_noeud(int num,FEM_NOEUD* noeud)
141     {
142     FEM_TEMPLATE_ELEMENT<3>::change_noeud(num,noeud);
143     }
144    
145     BOITE_3D& FEM_TRIANGLE3::get_boite_3D(void)
146     {
147     return FEM_TEMPLATE_ELEMENT<3>::get_boite_3D();
148     }
149    
150    
151     void FEM_TRIANGLE3::enregistrer(std::ostream& o)
152     {
153 francois 378 if (maillage!=NULL)
154     if (get_lien_topologie()!=NULL) o << "%" << get_id() << "=FEM_TRIANGLE3($"<< get_lien_topologie()->get_id() << ",$" << maillage->get_id() << ",$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<< ");" << std::endl;
155     else o << "%" << get_id() << "=FEM_TRIANGLE3(NULL,$" << maillage->get_id() << ",$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<< ");" << std::endl;
156     else
157     if (get_lien_topologie()!=NULL) o << "%" << get_id() << "=FEM_TRIANGLE3($"<< get_lien_topologie()->get_id() << ",NULL,$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<< ");" << std::endl;
158     else o << "%" << get_id() << "=FEM_TRIANGLE3(NULL,NULL,$" << tab[0]->get_id() << ",$" << tab[1]->get_id() << ",$" << tab[2]->get_id()<< ");" << std::endl;
159    
160 francois 637 }
161 francois 310
162 francois 637 int FEM_TRIANGLE3::get_nb_pt_gauss(int degre)
163     {
164     if (degre<2) return 1;
165     if (degre<3) return 3;
166     if (degre<4) return 4;
167     if (degre<5) return 6;
168     if (degre<6) return 7;
169     return 0;
170     }
171    
172     void FEM_TRIANGLE3::get_pt_gauss(int degre,int num,double &w,double *uv)
173     {
174     if (degre<2)
175     {
176     if (num==0) {w=0.5;uv[0]=0.333333333333333;uv[1]=0.333333333333333;}
177     return;
178     }
179     if (degre<3)
180     {
181     if (num==0) {w=0.166666666666667;uv[0]=0.166666666666667;uv[1]=0.166666666666667;}
182     if (num==1) {w=0.166666666666667;uv[0]=0.666666666666667;uv[1]=0.166666666666667;}
183     if (num==2) {w=0.166666666666667;uv[0]=0.166666666666667;uv[1]=0.666666666666667;}
184     return;
185     }
186     if (degre<4)
187     {
188     if (num==0) {w=-0.28125;uv[0]=0.333333333333333;uv[1]=0.333333333333333;}
189     if (num==1) {w=0.260416666666667;uv[0]=0.2;uv[1]=0.2;}
190     if (num==2) {w=0.260416666666667;uv[0]=0.6;uv[1]=0.2;}
191     if (num==3) {w=0.260416666666667;uv[0]=0.2;uv[1]=0.6;}
192     return;
193     }
194     if (degre<5)
195     {
196     if (num==0) {w=0.111690794839005;uv[0]=0.445948490915965;uv[1]=0.445948490915965;}
197     if (num==1) {w=0.111690794839005;uv[0]=0.108103018168070;uv[1]=0.445948490915965;}
198     if (num==2) {w=0.111690794839005;uv[0]=0.445948490915965;uv[1]=0.108103018168070;}
199     if (num==3) {w=0.054975871827661;uv[0]=0.091576213509771;uv[1]=0.091576213509771;}
200     if (num==4) {w=0.054975871827661;uv[0]=0.816847572980458;uv[1]=0.091576213509771;}
201     if (num==5) {w=0.054975871827661;uv[0]=0.091576213509771;uv[1]=0.816847572980458;}
202     return;
203     }
204     if (degre<6)
205     {
206     if (num==0) {w=0.1125;uv[0]=0.333333333333333;uv[1]=0.333333333333333;}
207     if (num==1) {w=0.066197076394253;uv[0]=0.470142064105115;uv[1]=0.470142064105115;}
208     if (num==2) {w=0.066197076394253;uv[0]=0.059715871789770;uv[1]=0.470142064105115;}
209     if (num==3) {w=0.066197076394253;uv[0]=0.470142064105115;uv[1]=0.059715871789770;}
210     if (num==4) {w=0.062969590272413;uv[0]=0.101286507323456;uv[1]=0.101286507323456;}
211     if (num==5) {w=0.062969590272413;uv[0]=0.797426985353088;uv[1]=0.101286507323456;}
212     if (num==6) {w=0.062969590272413;uv[0]=0.101286507323456;uv[1]=0.797426985353088;}
213     return;
214     }
215     }
216    
217 francois 638 int FEM_TRIANGLE3::get_nb_fonction_interpolation(void)
218 francois 310 {
219     return 3;
220     }
221 francois 638 int FEM_TRIANGLE3::get_degremax_fonction_interpolation(void)
222     {
223     return 1;
224     }
225 francois 310 double FEM_TRIANGLE3::get_fonction_interpolation(int num,double *uv)
226     {
227     double val;
228     switch (num)
229     {
230     case 1:
231     val=1-uv[0]-uv[1];
232     break;
233     case 2:
234     val=uv[0];
235     break;
236     case 3:
237     val=uv[1];
238     break;
239    
240     }
241     return val;
242     }
243    
244     double FEM_TRIANGLE3::get_fonction_derive_interpolation(int num,int num_variable,double *uv)
245     {
246     double val;
247    
248     switch (num)
249     {
250     case 1:
251     switch (num_variable)
252     {
253     case 1:
254     val=-1;
255     break;
256     case 2:
257     val=-1;
258     break;
259     } break;
260     case 2:
261     switch (num_variable)
262     {
263     case 1:
264     val=1;
265     break;
266     case 2:
267     val=0.;
268     break;
269     }break;
270     case 3:
271     switch (num_variable)
272     {
273     case 1:
274     val=0.;
275     break;
276     case 2:
277     val=1.;
278     break;
279     }break;
280    
281     }
282     return val;
283     }
284    
285     double FEM_TRIANGLE3::get_jacobien(double* jac,double *uv,int& li,int& col,double unite)
286 gervaislavoie 382 {
287 francois 647 int nb=get_nb_fem_noeud();
288     OT_VECTEUR_3D *point= new OT_VECTEUR_3D[nb];
289 gervaislavoie 382
290 francois 647 for (int i=0;i<nb;i++)
291     {
292     point[i].change_x(get_fem_noeud(i)->get_x());
293     point[i].change_y(get_fem_noeud(i)->get_y());
294     point[i].change_z(get_fem_noeud(i)->get_z());
295     }
296    
297     OT_MATRICE_3D J;
298     for (int i=0;i<2;i++)
299     {
300     J(i,0)=0.;J(i,1)=0.;J(i,2)=0.;
301     for (int k=0;k<nb;k++)
302     {
303     J(i,0)=J(i,0)+get_fonction_derive_interpolation(k+1,i+1,uv)*point[k].get_x();
304     J(i,1)=J(i,1)+get_fonction_derive_interpolation(k+1,i+1,uv)*point[k].get_y();
305     J(i,2)=J(i,2)+get_fonction_derive_interpolation(k+1,i+1,uv)*point[k].get_z();
306     }
307     }
308     OT_VECTEUR_3D J1(J(0,0),J(0,1),J(0,2));
309     OT_VECTEUR_3D J2(J(1,0),J(1,1),J(1,2));
310     OT_VECTEUR_3D J3=J1&J2;
311     J3.norme();
312     J(2,0)=J3(0);
313     J(2,1)=J3(1);
314     J(2,2)=J3(2);
315 gervaislavoie 389
316 francois 647 jac[0]=J(0,0);
317     jac[1]=J(0,1);
318     jac[2]=J(0,2);
319 gervaislavoie 389
320 francois 647 jac[3]=J(1,0);
321     jac[4]=J(1,1);
322     jac[5]=J(1,2);
323 gervaislavoie 389
324 francois 647 jac[6]=J(2,0);
325     jac[7]=J(2,1);
326     jac[8]=J(2,2);
327     li=3;col=3;
328 gervaislavoie 389
329 francois 647 return J.get_determinant();
330 francois 310 }
331    
332 francois 635
333     bool FEM_TRIANGLE3::valide_parametre_EF(double* uvw)
334     {
335     if (uvw[0]>=-1e-10)
336     if (uvw[1]>=-1e-10)
337     if (uvw[0]+uvw[1]<=1.+1e-10)
338     return true;
339     return false;
340     }
341    
342 francois 310 void FEM_TRIANGLE3::get_inverse_jacob(double* j,double *uv,double unite)
343     {
344 francois 647 double jac[9];
345     int li,col;
346     get_jacobien(jac,uv,li,col);
347     OT_MATRICE_3D J;
348     J(0,0)=jac[0];
349     J(0,1)=jac[1];
350     J(0,2)=jac[2];
351     J(1,0)=jac[3];
352     J(1,1)=jac[4];
353     J(1,2)=jac[5];
354     J(2,0)=jac[6];
355     J(2,1)=jac[7];
356     J(2,2)=jac[8];
357     OT_MATRICE_3D j_i=J.inverse();
358    
359 gervaislavoie 382
360     j[0]=j_i(0,0);
361 francois 647 j[1]=j_i(0,1);
362     j[2]=j_i(0,2);
363 gervaislavoie 382
364 francois 647 j[3]=j_i(1,0);
365 gervaislavoie 382 j[4]=j_i(1,1);
366 francois 647 j[5]=j_i(1,2);
367 gervaislavoie 382
368 francois 647 j[6]=j_i(2,0);
369     j[7]=j_i(2,1);
370 gervaislavoie 382 j[8]=j_i(2,2);
371 francois 310 }