ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/geometrie/src/fem_triangle3.cpp
Revision: 638
Committed: Mon Jan 26 21:56:20 2015 UTC (10 years, 3 months ago) by francois
File size: 16492 byte(s)
Log Message: 
ajout d'une méthode qui renvoit le degre max des fonctions d'interpolation et changement de nom de la fonction qui renvoie le nombre de fonction d'interpolation (ajout d'un get pour uniformiser)

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     // 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
288 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
289 gervaislavoie 382 FEM_NOEUD* n1=get_fem_noeud(0);
290     FEM_NOEUD* n2=get_fem_noeud(1);
291     FEM_NOEUD* n3=get_fem_noeud(2);
292    
293 gervaislavoie 389 double xyzn1[3];
294     double xyzn2[3];
295     double xyzn3[3];
296    
297 gervaislavoie 390 xyzn1[0]=n1->get_x();
298     xyzn1[1]=n1->get_y();
299     xyzn1[2]=n1->get_z();
300 gervaislavoie 389
301 gervaislavoie 390 xyzn2[0]=n2->get_x();
302     xyzn2[1]=n2->get_y();
303     xyzn2[2]=n2->get_z();
304 gervaislavoie 389
305 gervaislavoie 390 xyzn3[0]=n3->get_x();
306     xyzn3[1]=n3->get_y();
307     xyzn3[2]=n3->get_z();
308 gervaislavoie 389
309 gervaislavoie 392 // Position des noeuds dans le repère initial (i)
310     OT_VECTEUR_3D n1_i(xyzn1);
311     OT_VECTEUR_3D n2_i(xyzn2);
312     OT_VECTEUR_3D n3_i(xyzn3);
313     OT_MATRICE_3D pos_i(n1_i,n2_i,n3_i);
314    
315 gervaislavoie 382 // Le nouveau repère est défini par les vecteurs vec1 et vec2 et son origine est le noeud n1
316     OT_VECTEUR_3D vec1(xyzn1,xyzn2);
317     OT_VECTEUR_3D vec3(xyzn1,xyzn3);
318     vec1.norme();
319     OT_VECTEUR_3D n=vec1&vec3; // Normale au triangle au noeud n1
320     n.norme();
321     OT_VECTEUR_3D vec2=n&vec1; // Vecteur normal à vec1 dans le plan du triangle
322     vec2.norme();
323     OT_MATRICE_3D pif(vec1,vec2,n); // Matrice de passage du repère Ri au repère Rf
324    
325 gervaislavoie 392 // Position des noeuds dans le repère final (f)
326 gervaislavoie 382 OT_MATRICE_3D pif_inverse=pif.inverse();
327 gervaislavoie 392 OT_MATRICE_3D pos_f=pif_inverse*pos_i;
328     OT_VECTEUR_3D vecn1=pos_f.get_vecteur1();
329     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
330     pos_f.change_vecteur2(pos_f.get_vecteur2()-vecn1);
331     pos_f.change_vecteur3(pos_f.get_vecteur3()-vecn1);
332 gervaislavoie 382
333 gervaislavoie 389 // Jacobien dans le repère final
334     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));
335     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));
336     OT_VECTEUR_3D dNi_dZ(0.,0.,0.);
337 gervaislavoie 382
338 gervaislavoie 389 OT_MATRICE_3D dNi(dNi_de,dNi_dn,dNi_dZ);
339     OT_MATRICE_3D dNi_transpose=dNi.transpose();
340 gervaislavoie 392 OT_MATRICE_3D pos_f_transpose=pos_f.transpose();
341     OT_MATRICE_3D J_f=dNi_transpose*pos_f_transpose;
342 gervaislavoie 382
343 gervaislavoie 389 // Jacobien dans le repère initial
344     OT_MATRICE_3D J_i=pif*J_f;
345 gervaislavoie 382
346     jac[0]=J_i(0,0);
347     jac[1]=J_i(1,0);
348     jac[2]=J_i(2,0);
349    
350     jac[3]=J_i(0,1);
351     jac[4]=J_i(1,1);
352     jac[5]=J_i(2,1);
353    
354     jac[6]=0.;
355     jac[7]=0.;
356     jac[8]=0.;
357 francois 310 }
358    
359 francois 635
360     bool FEM_TRIANGLE3::valide_parametre_EF(double* uvw)
361     {
362     if (uvw[0]>=-1e-10)
363     if (uvw[1]>=-1e-10)
364     if (uvw[0]+uvw[1]<=1.+1e-10)
365     return true;
366     return false;
367     }
368    
369 francois 310 void FEM_TRIANGLE3::get_inverse_jacob(double* j,double *uv,double unite)
370     {
371 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
372 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
373 gervaislavoie 382 FEM_NOEUD* n1=get_fem_noeud(0);
374     FEM_NOEUD* n2=get_fem_noeud(1);
375     FEM_NOEUD* n3=get_fem_noeud(2);
376    
377 gervaislavoie 389 double xyzn1[3];
378     double xyzn2[3];
379     double xyzn3[3];
380 gervaislavoie 382
381 gervaislavoie 390 xyzn1[0]=n1->get_x();
382     xyzn1[1]=n1->get_y();
383     xyzn1[2]=n1->get_z();
384 gervaislavoie 382
385 gervaislavoie 390 xyzn2[0]=n2->get_x();
386     xyzn2[1]=n2->get_y();
387     xyzn2[2]=n2->get_z();
388 gervaislavoie 382
389 gervaislavoie 390 xyzn3[0]=n3->get_x();
390     xyzn3[1]=n3->get_y();
391     xyzn3[2]=n3->get_z();
392 gervaislavoie 392
393     // Position des noeuds dans le repère initial (i)
394     OT_VECTEUR_3D n1_i(xyzn1);
395     OT_VECTEUR_3D n2_i(xyzn2);
396     OT_VECTEUR_3D n3_i(xyzn3);
397     OT_MATRICE_3D pos_i(n1_i,n2_i,n3_i);
398 gervaislavoie 382
399 gervaislavoie 392 // Le nouveau repère est défini par les vecteurs vec1 et vec2 et son origine est le noeud n1
400 gervaislavoie 382 OT_VECTEUR_3D vec1(xyzn1,xyzn2);
401     OT_VECTEUR_3D vec3(xyzn1,xyzn3);
402     vec1.norme();
403     OT_VECTEUR_3D n=vec1&vec3; // Normale au triangle au noeud n1
404     n.norme();
405     OT_VECTEUR_3D vec2=n&vec1; // Vecteur normal à vec1 dans le plan du triangle
406     vec2.norme();
407     OT_MATRICE_3D pif(vec1,vec2,n); // Matrice de passage du repère Ri au repère Rf
408    
409 gervaislavoie 392 // Position des noeuds dans le repère final (f)
410 gervaislavoie 382 OT_MATRICE_3D pif_inverse=pif.inverse();
411 gervaislavoie 392 OT_MATRICE_3D pos_f=pif_inverse*pos_i;
412     OT_VECTEUR_3D vecn1=pos_f.get_vecteur1();
413     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
414     pos_f.change_vecteur2(pos_f.get_vecteur2()-vecn1);
415     pos_f.change_vecteur3(pos_f.get_vecteur3()-vecn1);
416 gervaislavoie 382
417 gervaislavoie 389 // Déterminant du jacobien 2D dans le repère final
418     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));
419     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));
420     OT_VECTEUR_3D dNi_dZ(0.,0.,0.);
421 gervaislavoie 382
422 gervaislavoie 389 OT_MATRICE_3D dNi(dNi_de,dNi_dn,dNi_dZ);
423     OT_MATRICE_3D dNi_transpose=dNi.transpose();
424 gervaislavoie 392 OT_MATRICE_3D pos_f_transpose=pos_f.transpose();
425     OT_MATRICE_3D J=dNi_transpose*pos_f_transpose;
426 gervaislavoie 382
427 gervaislavoie 389 double J11=J(0,0);
428     double J21=J(1,0);
429     double J12=J(0,1);
430     double J22=J(1,1);
431 gervaislavoie 382 double det_J=J11*J22-J12*J21;
432    
433 gervaislavoie 389 // Jacobien inverse 3D dans le repère final
434 gervaislavoie 382 OT_VECTEUR_3D col1_jf_3d(J22/det_J,-J21/det_J,0.);
435     OT_VECTEUR_3D col2_jf_3d(-J12/det_J,J11/det_J,0.);
436     OT_VECTEUR_3D col3_jf_3d(0.,0.,0.);
437 gervaislavoie 389 OT_MATRICE_3D j_f_3d(col1_jf_3d,col2_jf_3d,col3_jf_3d);
438 gervaislavoie 382
439 gervaislavoie 389 // Jacobien inverse 3D dans le repère initial
440     OT_MATRICE_3D j_i=pif*j_f_3d;
441 gervaislavoie 382
442     j[0]=j_i(0,0);
443     j[1]=j_i(1,0);
444     j[2]=j_i(2,0);
445    
446     j[3]=j_i(0,1);
447     j[4]=j_i(1,1);
448     j[5]=j_i(2,1);
449    
450     j[6]=j_i(0,2);
451     j[7]=j_i(1,2);
452     j[8]=j_i(2,2);
453 francois 310 }