ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_triangle3.cpp
Revision: 637
Committed: Mon Jan 26 21:26:37 2015 UTC (10 years, 3 months ago) by francois
File size: 16407 byte(s)
Log Message:
entree des points de gauss pour les différents éléments finis depuis reference JCC

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 310 int FEM_TRIANGLE3::nb_fonction_interpolation(void)
218     {
219     return 3;
220     }
221    
222     double FEM_TRIANGLE3::get_fonction_interpolation(int num,double *uv)
223     {
224     double val;
225     switch (num)
226     {
227     case 1:
228     val=1-uv[0]-uv[1];
229     break;
230     case 2:
231     val=uv[0];
232     break;
233     case 3:
234     val=uv[1];
235     break;
236    
237     }
238     return val;
239     }
240    
241     double FEM_TRIANGLE3::get_fonction_derive_interpolation(int num,int num_variable,double *uv)
242     {
243     double val;
244    
245     switch (num)
246     {
247     case 1:
248     switch (num_variable)
249     {
250     case 1:
251     val=-1;
252     break;
253     case 2:
254     val=-1;
255     break;
256     } break;
257     case 2:
258     switch (num_variable)
259     {
260     case 1:
261     val=1;
262     break;
263     case 2:
264     val=0.;
265     break;
266     }break;
267     case 3:
268     switch (num_variable)
269     {
270     case 1:
271     val=0.;
272     break;
273     case 2:
274     val=1.;
275     break;
276     }break;
277    
278     }
279     return val;
280     }
281    
282     double FEM_TRIANGLE3::get_jacobien(double* jac,double *uv,int& li,int& col,double unite)
283 gervaislavoie 382 {
284     // Il faut exprimer les coordonnées x,y,z du repère 3D initial Ri dans un repère 2D final Rf de coordonnées
285 gervaislavoie 389 // x',y',z' où z' est nul afin de pouvoir utiliser les fonctions d'interpolation de l'élément de référence 2D
286 gervaislavoie 382 FEM_NOEUD* n1=get_fem_noeud(0);
287     FEM_NOEUD* n2=get_fem_noeud(1);
288     FEM_NOEUD* n3=get_fem_noeud(2);
289    
290 gervaislavoie 389 double xyzn1[3];
291     double xyzn2[3];
292     double xyzn3[3];
293    
294 gervaislavoie 390 xyzn1[0]=n1->get_x();
295     xyzn1[1]=n1->get_y();
296     xyzn1[2]=n1->get_z();
297 gervaislavoie 389
298 gervaislavoie 390 xyzn2[0]=n2->get_x();
299     xyzn2[1]=n2->get_y();
300     xyzn2[2]=n2->get_z();
301 gervaislavoie 389
302 gervaislavoie 390 xyzn3[0]=n3->get_x();
303     xyzn3[1]=n3->get_y();
304     xyzn3[2]=n3->get_z();
305 gervaislavoie 389
306 gervaislavoie 392 // Position des noeuds dans le repère initial (i)
307     OT_VECTEUR_3D n1_i(xyzn1);
308     OT_VECTEUR_3D n2_i(xyzn2);
309     OT_VECTEUR_3D n3_i(xyzn3);
310     OT_MATRICE_3D pos_i(n1_i,n2_i,n3_i);
311    
312 gervaislavoie 382 // Le nouveau repère est défini par les vecteurs vec1 et vec2 et son origine est le noeud n1
313     OT_VECTEUR_3D vec1(xyzn1,xyzn2);
314     OT_VECTEUR_3D vec3(xyzn1,xyzn3);
315     vec1.norme();
316     OT_VECTEUR_3D n=vec1&vec3; // Normale au triangle au noeud n1
317     n.norme();
318     OT_VECTEUR_3D vec2=n&vec1; // Vecteur normal à vec1 dans le plan du triangle
319     vec2.norme();
320     OT_MATRICE_3D pif(vec1,vec2,n); // Matrice de passage du repère Ri au repère Rf
321    
322 gervaislavoie 392 // Position des noeuds dans le repère final (f)
323 gervaislavoie 382 OT_MATRICE_3D pif_inverse=pif.inverse();
324 gervaislavoie 392 OT_MATRICE_3D pos_f=pif_inverse*pos_i;
325     OT_VECTEUR_3D vecn1=pos_f.get_vecteur1();
326     pos_f.change_vecteur1(pos_f.get_vecteur1()-vecn1); // Étant donné que le noeud 1 est l'origine du repère final, on soustrait ses coordonnées à chaque noeud
327     pos_f.change_vecteur2(pos_f.get_vecteur2()-vecn1);
328     pos_f.change_vecteur3(pos_f.get_vecteur3()-vecn1);
329 gervaislavoie 382
330 gervaislavoie 389 // Jacobien dans le repère final
331     OT_VECTEUR_3D dNi_de(get_fonction_derive_interpolation(1,1,uv),get_fonction_derive_interpolation(2,1,uv),get_fonction_derive_interpolation(3,1,uv));
332     OT_VECTEUR_3D dNi_dn(get_fonction_derive_interpolation(1,2,uv),get_fonction_derive_interpolation(2,2,uv),get_fonction_derive_interpolation(3,2,uv));
333     OT_VECTEUR_3D dNi_dZ(0.,0.,0.);
334 gervaislavoie 382
335 gervaislavoie 389 OT_MATRICE_3D dNi(dNi_de,dNi_dn,dNi_dZ);
336     OT_MATRICE_3D dNi_transpose=dNi.transpose();
337 gervaislavoie 392 OT_MATRICE_3D pos_f_transpose=pos_f.transpose();
338     OT_MATRICE_3D J_f=dNi_transpose*pos_f_transpose;
339 gervaislavoie 382
340 gervaislavoie 389 // Jacobien dans le repère initial
341     OT_MATRICE_3D J_i=pif*J_f;
342 gervaislavoie 382
343     jac[0]=J_i(0,0);
344     jac[1]=J_i(1,0);
345     jac[2]=J_i(2,0);
346    
347     jac[3]=J_i(0,1);
348     jac[4]=J_i(1,1);
349     jac[5]=J_i(2,1);
350    
351     jac[6]=0.;
352     jac[7]=0.;
353     jac[8]=0.;
354 francois 310 }
355    
356 francois 635
357     bool FEM_TRIANGLE3::valide_parametre_EF(double* uvw)
358     {
359     if (uvw[0]>=-1e-10)
360     if (uvw[1]>=-1e-10)
361     if (uvw[0]+uvw[1]<=1.+1e-10)
362     return true;
363     return false;
364     }
365    
366 francois 310 void FEM_TRIANGLE3::get_inverse_jacob(double* j,double *uv,double unite)
367     {
368 gervaislavoie 382 // Il faut exprimer les coordonnées x,y,z du repère 3D initial Ri dans un repère 2D final Rf de coordonnées
369 gervaislavoie 389 // x',y',z' où z' est nul afin de pouvoir utiliser les fonctions d'interpolation de l'élément de référence
370 gervaislavoie 382 FEM_NOEUD* n1=get_fem_noeud(0);
371     FEM_NOEUD* n2=get_fem_noeud(1);
372     FEM_NOEUD* n3=get_fem_noeud(2);
373    
374 gervaislavoie 389 double xyzn1[3];
375     double xyzn2[3];
376     double xyzn3[3];
377 gervaislavoie 382
378 gervaislavoie 390 xyzn1[0]=n1->get_x();
379     xyzn1[1]=n1->get_y();
380     xyzn1[2]=n1->get_z();
381 gervaislavoie 382
382 gervaislavoie 390 xyzn2[0]=n2->get_x();
383     xyzn2[1]=n2->get_y();
384     xyzn2[2]=n2->get_z();
385 gervaislavoie 382
386 gervaislavoie 390 xyzn3[0]=n3->get_x();
387     xyzn3[1]=n3->get_y();
388     xyzn3[2]=n3->get_z();
389 gervaislavoie 392
390     // Position des noeuds dans le repère initial (i)
391     OT_VECTEUR_3D n1_i(xyzn1);
392     OT_VECTEUR_3D n2_i(xyzn2);
393     OT_VECTEUR_3D n3_i(xyzn3);
394     OT_MATRICE_3D pos_i(n1_i,n2_i,n3_i);
395 gervaislavoie 382
396 gervaislavoie 392 // Le nouveau repère est défini par les vecteurs vec1 et vec2 et son origine est le noeud n1
397 gervaislavoie 382 OT_VECTEUR_3D vec1(xyzn1,xyzn2);
398     OT_VECTEUR_3D vec3(xyzn1,xyzn3);
399     vec1.norme();
400     OT_VECTEUR_3D n=vec1&vec3; // Normale au triangle au noeud n1
401     n.norme();
402     OT_VECTEUR_3D vec2=n&vec1; // Vecteur normal à vec1 dans le plan du triangle
403     vec2.norme();
404     OT_MATRICE_3D pif(vec1,vec2,n); // Matrice de passage du repère Ri au repère Rf
405    
406 gervaislavoie 392 // Position des noeuds dans le repère final (f)
407 gervaislavoie 382 OT_MATRICE_3D pif_inverse=pif.inverse();
408 gervaislavoie 392 OT_MATRICE_3D pos_f=pif_inverse*pos_i;
409     OT_VECTEUR_3D vecn1=pos_f.get_vecteur1();
410     pos_f.change_vecteur1(pos_f.get_vecteur1()-vecn1); // Étant donné que le noeud 1 est l'origine du repère final, on soustrait ses coordonnées à chaque noeud
411     pos_f.change_vecteur2(pos_f.get_vecteur2()-vecn1);
412     pos_f.change_vecteur3(pos_f.get_vecteur3()-vecn1);
413 gervaislavoie 382
414 gervaislavoie 389 // Déterminant du jacobien 2D dans le repère final
415     OT_VECTEUR_3D dNi_de(get_fonction_derive_interpolation(1,1,uv),get_fonction_derive_interpolation(2,1,uv),get_fonction_derive_interpolation(3,1,uv));
416     OT_VECTEUR_3D dNi_dn(get_fonction_derive_interpolation(1,2,uv),get_fonction_derive_interpolation(2,2,uv),get_fonction_derive_interpolation(3,2,uv));
417     OT_VECTEUR_3D dNi_dZ(0.,0.,0.);
418 gervaislavoie 382
419 gervaislavoie 389 OT_MATRICE_3D dNi(dNi_de,dNi_dn,dNi_dZ);
420     OT_MATRICE_3D dNi_transpose=dNi.transpose();
421 gervaislavoie 392 OT_MATRICE_3D pos_f_transpose=pos_f.transpose();
422     OT_MATRICE_3D J=dNi_transpose*pos_f_transpose;
423 gervaislavoie 382
424 gervaislavoie 389 double J11=J(0,0);
425     double J21=J(1,0);
426     double J12=J(0,1);
427     double J22=J(1,1);
428 gervaislavoie 382 double det_J=J11*J22-J12*J21;
429    
430 gervaislavoie 389 // Jacobien inverse 3D dans le repère final
431 gervaislavoie 382 OT_VECTEUR_3D col1_jf_3d(J22/det_J,-J21/det_J,0.);
432     OT_VECTEUR_3D col2_jf_3d(-J12/det_J,J11/det_J,0.);
433     OT_VECTEUR_3D col3_jf_3d(0.,0.,0.);
434 gervaislavoie 389 OT_MATRICE_3D j_f_3d(col1_jf_3d,col2_jf_3d,col3_jf_3d);
435 gervaislavoie 382
436 gervaislavoie 389 // Jacobien inverse 3D dans le repère initial
437     OT_MATRICE_3D j_i=pif*j_f_3d;
438 gervaislavoie 382
439     j[0]=j_i(0,0);
440     j[1]=j_i(1,0);
441     j[2]=j_i(2,0);
442    
443     j[3]=j_i(0,1);
444     j[4]=j_i(1,1);
445     j[5]=j_i(2,1);
446    
447     j[6]=j_i(0,2);
448     j[7]=j_i(1,2);
449     j[8]=j_i(2,2);
450 francois 310 }