1 #ifndef MEDMEM_INTERPOLATION_TOOLS_HXX
2 #define MEDMEM_INTERPOLATION_TOOLS_HXX
4 #include "MEDMEM_define.hxx"
6 #define _TEMPLATE_ template <class CHAMP,class VALEURCHAMP,class NUAGENOEUD, class NOEUD, class NUAGEMAILLE>
7 #define _PARAM_ CHAMP,VALEURCHAMP,NUAGENOEUD,NOEUD,NUAGEMAILLE
9 #define face2367(x,y,z) ((x6*(-y2 + y3) + x3*(y2 - y6) + x2*(-y3 + y6))*z - x6*y3*z2 + x3*y6*z2 + x6*y2*z3 - x2*y6*z3 - x3*y2*z6 + x2*y3*z6 + y*(x6*(z2 - z3) + x2*(z3 - z6) + x3*(-z2 + z6)) + x*(y6*(-z2 + z3) + y3*(z2 - z6) + y2*(-z3 + z6)))
10 #define face4567(x,y,z) ((x6*(-y4 + y5) + x5*(y4 - y6) + x4*(-y5 + y6))*z - x6*y5*z4 + x5*y6*z4 + x6*y4*z5 - x4*y6*z5 - x5*y4*z6 + x4*y5*z6 + y*(x6*(z4 - z5) + x4*(z5 - z6) + x5*(-z4 + z6)) + x*(y6*(-z4 + z5) + y5*(z4 - z6) + y4*(-z5 + z6)))
11 #define face1256(x,y,z) ((x5*(-y1 + y2) + x2*(y1 - y5) + x1*(-y2 + y5))*z - x5*y2*z1 + x2*y5*z1 + x5*y1*z2 - x1*y5*z2 - x2*y1*z5 + x1*y2*z5 + y*(x5*(z1 - z2) + x1*(z2 - z5) + x2*(-z1 + z5)) + x*(y5*(-z1 + z2) + y2*(z1 - z5) + y1*(-z2 + z5)))
12 #define face0347(x,y,z) ((x4*(-y0 + y3) + x3*(y0 - y4) + x0*(-y3 + y4))*z - x4*y3*z0 + x3*y4*z0 + x4*y0*z3 - x0*y4*z3 - x3*y0*z4 + x0*y3*z4 + y*(x4*(z0 - z3) + x0*(z3 - z4) + x3*(-z0 + z4)) + x*(y4*(-z0 + z3) + y3*(z0 - z4) + y0*(-z3 + z4)))
13 #define face0145(x,y,z) ((x4*(-y0 + y1) + x1*(y0 - y4) + x0*(-y1 + y4))*z - x4*y1*z0 + x1*y4*z0 + x4*y0*z1 - x0*y4*z1 - x1*y0*z4 + x0*y1*z4 + y*(x4*(z0 - z1) + x0*(z1 - z4) + x1*(-z0 + z4)) + x*(y4*(-z0 + z1) + y1*(z0 - z4) + y0*(-z1 + z4)))
14 #define face0123(x,y,z) ((x2*(-y0 + y1) + x1*(y0 - y2) + x0*(-y1 + y2))*z - x2*y1*z0 + x1*y2*z0 + x2*y0*z1 - x0*y2*z1 - x1*y0*z2 + x0*y1*z2 + y*(x2*(z0 - z1) + x0*(z1 - z2) + x1*(-z0 + z2)) + x*(y2*(-z0 + z1) + y1*(z0 - z2) + y0*(-z1 + z2)))
20 _TEMPLATE_ class Calcul_Interpolation
24 NUAGEMAILLE * mailles;
27 Calcul_Interpolation():noeuds(NULL),mailles(NULL),champ(NULL) {}
28 Calcul_Interpolation(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):noeuds(nn),mailles(nm),champ(c) {}
29 ~Calcul_Interpolation() {}
30 virtual VALEURCHAMP operator() (const NOEUD & n, int num_maille) {cerr<<"APPEL OPERATOR() DE LA CLASSE MERE CALCUL_INTERPOLATION => EXIT(-1)"<<endl;exit(-1);}
33 _TEMPLATE_ class Calcul_Interpolation_P0;
35 _TEMPLATE_ class Calcul_Interpolation_Tria3;
36 _TEMPLATE_ class Calcul_Interpolation_Tria6;
37 _TEMPLATE_ class Calcul_Interpolation_Quad4;
38 _TEMPLATE_ class Calcul_Interpolation_Quad8;
39 _TEMPLATE_ class Calcul_Interpolation_Tetra4;
40 _TEMPLATE_ class Calcul_Interpolation_Tetra10;
41 _TEMPLATE_ class Calcul_Interpolation_Hexa8;
42 _TEMPLATE_ class Calcul_Interpolation_Hexa20;
43 _TEMPLATE_ class Calcul_Interpolation_Penta6;
44 _TEMPLATE_ class Calcul_Interpolation_Penta15;
45 _TEMPLATE_ class Calcul_Interpolation_Pyra5;
46 _TEMPLATE_ class Calcul_Interpolation_Pyra13;
48 _TEMPLATE_ class Calcul_Hybride
51 NUAGEMAILLE * mailles;
52 map<int,Calcul_Interpolation<_PARAM_> *> fonctions;
54 Calcul_Hybride():mailles(NULL) {}
55 Calcul_Hybride(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c);
57 VALEURCHAMP operator() (const NOEUD & n, int num_maille);
62 _TEMPLATE_ Calcul_Hybride<_PARAM_>::Calcul_Hybride(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):mailles(nm)
64 fonctions[MED_TRIA3 ]=new Calcul_Interpolation_Tria3 <_PARAM_>(nn,nm,c);
65 fonctions[MED_TRIA6 ]=new Calcul_Interpolation_Tria6 <_PARAM_>(nn,nm,c);
66 fonctions[MED_QUAD4 ]=new Calcul_Interpolation_Quad4 <_PARAM_>(nn,nm,c);
67 fonctions[MED_QUAD8 ]=new Calcul_Interpolation_Quad8 <_PARAM_>(nn,nm,c);
68 fonctions[MED_TETRA4 ]=new Calcul_Interpolation_Tetra4 <_PARAM_>(nn,nm,c);
69 fonctions[MED_TETRA10]=new Calcul_Interpolation_Tetra10<_PARAM_>(nn,nm,c);
70 fonctions[MED_HEXA8 ]=new Calcul_Interpolation_Hexa8 <_PARAM_>(nn,nm,c);
71 fonctions[MED_HEXA20 ]=new Calcul_Interpolation_Hexa20 <_PARAM_>(nn,nm,c);
72 fonctions[MED_PENTA6 ]=new Calcul_Interpolation_Penta6 <_PARAM_>(nn,nm,c);
73 fonctions[MED_PENTA15]=new Calcul_Interpolation_Penta15<_PARAM_>(nn,nm,c);
74 fonctions[MED_PYRA5 ]=new Calcul_Interpolation_Pyra5 <_PARAM_>(nn,nm,c);
75 fonctions[MED_PYRA13 ]=new Calcul_Interpolation_Pyra13 <_PARAM_>(nn,nm,c);
78 _TEMPLATE_ VALEURCHAMP Calcul_Hybride<_PARAM_>::operator() (const NOEUD & n, int num_maille)
80 return fonctions[mailles->DONNE_TYPE_MAILLE(num_maille)]->operator()(n,num_maille);
83 _TEMPLATE_ class Calcul_Interpolation_P0 : public Calcul_Interpolation<_PARAM_>
85 public : Calcul_Interpolation_P0(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
86 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
88 return (*champ)[num_maille];
91 _TEMPLATE_ class Calcul_Interpolation_Tria3 : public Calcul_Interpolation<_PARAM_>
93 public : Calcul_Interpolation_Tria3(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
94 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
96 int num0=mailles->DONNE_SOMMET_MAILLE(num_maille,0);
97 int num1=mailles->DONNE_SOMMET_MAILLE(num_maille,1);
98 int num2=mailles->DONNE_SOMMET_MAILLE(num_maille,2);
100 double x0=(*noeuds)[num0][0];
101 double y0=(*noeuds)[num0][1];
102 double x1=(*noeuds)[num1][0];
103 double y1=(*noeuds)[num1][1];
104 double x2=(*noeuds)[num2][0];
105 double y2=(*noeuds)[num2][1];
107 VALEURCHAMP v0=(*champ)[num0];
108 VALEURCHAMP v1=(*champ)[num1];
109 VALEURCHAMP v2=(*champ)[num2];
114 double lambda0=(y1-y2)*x+(x2-x1)*y+(x1*y2-x2*y1);
115 double lambda1=(y2-y0)*x+(x0-x2)*y+(x2*y0-x0*y2);
116 double lambda2=(y0-y1)*x+(x1-x0)*y+(x0*y1-x1*y0);
118 double delta = (x2-x1)*y0+(x0-x2)*y1+(x1-x0)*y2;
120 VALEURCHAMP retour(v0.SIZE());
122 return (1/delta)*(lambda0*v0+lambda1*v1+lambda2*v2);
127 _TEMPLATE_ class Calcul_Interpolation_Tria6 : public Calcul_Interpolation<_PARAM_>
129 public : Calcul_Interpolation_Tria6(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
130 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
132 // ON SUPPOSE IMPLICITEMENT QUE LES NOEUDS SUPPLEMENTAIRES SONT BIEN DES NOEUDS MILIEUX
133 int num0=mailles->DONNE_SOMMET_MAILLE(num_maille,0);
134 int num1=mailles->DONNE_SOMMET_MAILLE(num_maille,1);
135 int num2=mailles->DONNE_SOMMET_MAILLE(num_maille,2);
136 int num01=mailles->DONNE_SOMMET_MAILLE(num_maille,3);
137 int num12=mailles->DONNE_SOMMET_MAILLE(num_maille,4);
138 int num20=mailles->DONNE_SOMMET_MAILLE(num_maille,5);
140 double x0=(*noeuds)[num0][0];
141 double y0=(*noeuds)[num0][1];
142 double x1=(*noeuds)[num1][0];
143 double y1=(*noeuds)[num1][1];
144 double x2=(*noeuds)[num2][0];
145 double y2=(*noeuds)[num2][1];
147 VALEURCHAMP v0=(*champ)[num0];
148 VALEURCHAMP v1=(*champ)[num1];
149 VALEURCHAMP v2=(*champ)[num2];
150 VALEURCHAMP v01=(*champ)[num01];
151 VALEURCHAMP v12=(*champ)[num12];
152 VALEURCHAMP v20=(*champ)[num20];
157 double lambda0=(y1-y2)*x+(x2-x1)*y+(x1*y2-x2*y1);
158 double lambda1=(y2-y0)*x+(x0-x2)*y+(x2*y0-x0*y2);
159 double lambda2=(y0-y1)*x+(x1-x0)*y+(x0*y1-x1*y0);
161 double delta = (x2-x1)*y0+(x0-x2)*y1+(x1-x0)*y2;
164 // VALEURCHAMP retour(v0.SIZE()); //
166 return 2*(lambda0*lambda0*v0+
169 2*(lambda0*lambda1*v01+
171 lambda2*lambda0*v20))/(delta*delta)+
172 (lambda0*v0+lambda1*v1+lambda2*v2)/(-delta);
177 _TEMPLATE_ class Calcul_Interpolation_Quad4 : public Calcul_Interpolation<_PARAM_>
179 public : Calcul_Interpolation_Quad4(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
180 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
182 int num0=mailles->DONNE_SOMMET_MAILLE(num_maille,0);
183 int num1=mailles->DONNE_SOMMET_MAILLE(num_maille,1);
184 int num2=mailles->DONNE_SOMMET_MAILLE(num_maille,2);
185 int num3=mailles->DONNE_SOMMET_MAILLE(num_maille,3);
187 double x0=(*noeuds)[num0][0];
188 double y0=(*noeuds)[num0][1];
189 double x1=(*noeuds)[num1][0];
190 double y1=(*noeuds)[num1][1];
191 double x2=(*noeuds)[num2][0];
192 double y2=(*noeuds)[num2][1];
193 double x3=(*noeuds)[num3][0];
194 double y3=(*noeuds)[num3][1];
196 VALEURCHAMP v0=(*champ)[num0];
197 VALEURCHAMP v1=(*champ)[num1];
198 VALEURCHAMP v2=(*champ)[num2];
199 VALEURCHAMP v3=(*champ)[num3];
205 double mu0=-((x3*(y1-y2)+x1*(y2-y3)+x2*(y3-y1))*x*y+(x2*y2*(y1-y3)+x3*(-y1+y2)*y3+x1*y1*(y3-y2))*x+(x2*x3*(y2-y3)+x1*(x2*(y1-y2)+x3*(y3-y1)))*y+(x2*x3*y1*(y3-y2)+x1*(x3*y2*(y1-y3)+x2*(y2-y1)*y3)));
206 double mu1=(x0*(y2-y3)+x2*(y3-y0)+x3*(y0-y2))*x*y+(x3*y3*(y2-y0)+x0*(-y2+y3)*y0+x2*y2*(y0-y3))*x+(x3*x0*(y3-y0)+x2*(x3*(y2-y3)+x0*(y0-y2)))*y+(x2*x0*y2*(y0-y2)+x2*(x0*y2*(y2-y0)+x3*(y3-y2)*y0));
207 double mu2=-((x1*(y3-y0)+x3*(y0-y1)+x0*(y1-y3))*x*y+(x0*y0*(y3-y1)+x1*(-y3+y0)*y1+x3*y3*(y1-y0))*x+(x0*x1*(y0-y1)+x3*(x0*(y3-y0)+x1*(y1-y3)))*y+(x3*x1*y3*(y1-y2)+x3*(x1*y2*(y3-y1)+x0*(y0-y3)*y1)));
208 double mu3=(x2*(y0-y1)+x0*(y1-y2)+x1*(y2-y0))*x*y+(x1*y1*(y0-y2)+x2*(-y0+y1)*y2+x0*y0*(y2-y1))*x+(x1*x2*(y1-y2)+x0*(x1*(y0-y1)+x2*(y2-y0)))*y+(x0*x2*y0*(y2-y2)+x0*(x2*y2*(y0-y2)+x1*(y1-y0)*y2));
210 double delta=(y0-y2)*(y1-y3)*(x0*x2+x1*x3)-(x1*x2+x0*x3)*(y1-y2)*(y0-y3)-(x0*x1+x2*x3)*(y0-y1)*(y2-y3);
213 cout<<" ### Pour ( "<<x<<" , "<<y<<" )"<<endl;
214 cout<<" ### delta = "<<delta<<endl;
215 cout<<" ### Mu0 = "<<mu0<<endl;
216 cout<<" ### Mu1 = "<<mu1<<endl;
217 cout<<" ### Mu2 = "<<mu2<<endl;
218 cout<<" ### Mu3 = "<<mu3<<endl;
220 // VALEURCHAMP retour(v0.SIZE()); //
222 return (mu0*v0+mu1*v1+mu2*v2+mu3*v3)/delta;
227 _TEMPLATE_ class Calcul_Interpolation_Quad8 : public Calcul_Interpolation<_PARAM_>
229 public : Calcul_Interpolation_Quad8(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
230 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
232 cerr<<"Interpolation Q2 pas encore implémentée"<<endl;
236 _TEMPLATE_ class Calcul_Interpolation_Tetra4 : public Calcul_Interpolation<_PARAM_>
238 public : Calcul_Interpolation_Tetra4(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
239 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
241 int num0=mailles->DONNE_SOMMET_MAILLE(num_maille,0);
242 int num1=mailles->DONNE_SOMMET_MAILLE(num_maille,1);
243 int num2=mailles->DONNE_SOMMET_MAILLE(num_maille,2);
244 int num3=mailles->DONNE_SOMMET_MAILLE(num_maille,3);
246 double x0=(*noeuds)[num0][0];
247 double y0=(*noeuds)[num0][1];
248 double z0=(*noeuds)[num0][2];
249 double x1=(*noeuds)[num1][0];
250 double y1=(*noeuds)[num1][1];
251 double z1=(*noeuds)[num1][2];
252 double x2=(*noeuds)[num2][0];
253 double y2=(*noeuds)[num2][1];
254 double z2=(*noeuds)[num2][2];
255 double x3=(*noeuds)[num3][0];
256 double y3=(*noeuds)[num3][1];
257 double z3=(*noeuds)[num3][2];
259 VALEURCHAMP v0=(*champ)[num0];
260 VALEURCHAMP v1=(*champ)[num1];
261 VALEURCHAMP v2=(*champ)[num2];
262 VALEURCHAMP v3=(*champ)[num3];
268 double lambda0=(y3*(z1-z2)+y1*(z2-z3)+y2*(z3-z1))*x+(x3*(z2-z1)+x2*(z1-z3)+x1*(z3-z2))*y+(x3*(y1-y2)+x1*(y2-y3)+x2*(y3-y1))*z+((x3*y2-x2*y3)*z1+(x1*y3-x3*y1)*z2+(x2*y1-x1*y2)*z3);
269 double lambda1=-((y0*(z2-z3)+y2*(z3-z0)+y3*(z0-z2))*x+(x0*(z3-z2)+x3*(z2-z0)+x2*(z0-z3))*y+(x0*(y2-y3)+x2*(y3-y0)+x3*(y0-y2))*z+((x0*y3-x2*y0)*z2+(x2*y0-x0*y2)*z3+(x3*y2-x2*y3)*z0));
270 double lambda2=(y1*(z3-z0)+y3*(z0-z1)+y0*(z1-z3))*x+(x1*(z0-z3)+x0*(z3-z1)+x3*(z1-z0))*y+(x1*(y3-y0)+x3*(y0-y1)+x0*(y1-y3))*z+((x1*y0-x2*y1)*z3+(x3*y1-x1*y3)*z0+(x0*y3-x3*y0)*z1);
271 double lambda3=-((y2*(z0-z1)+y0*(z1-z2)+y1*(z2-z0))*x+(x2*(z1-z0)+x1*(z0-z2)+x0*(z2-z1))*y+(x2*(y0-y1)+x0*(y1-y2)+x1*(y2-y0))*z+((x2*y1-x2*y2)*z0+(x0*y2-x2*y0)*z1+(x1*y0-x0*y1)*z2));
273 // A SIMPLFIFIER, FEIGNASSE !
274 double delta = (x1*y2-x1*y3)*z0 +(x0-y3- x0*y2)*z1 + x0*y3*z1 - x1*y0*z2 + x0*y1*z2 - x0*y3*z2 + x1*y3*z2 + x3*(y1*z0 - y2*z0 - y0*z1 + y2*z1 + y0*z2 - y1*z2) + x1*y0*z3 - x0*y1*z3 + x0*y2*z3 - x1*y2*z3 + x2*(y3*z0 + y0*z1 - y3*z1 - y0*z3 + y1*(-z0 + z3));
276 VALEURCHAMP retour(v0.SIZE()); //
278 retour=(lambda0*v0+lambda1*v1+lambda2*v2+lambda3*v3)/delta;
283 _TEMPLATE_ class Calcul_Interpolation_Tetra10 : public Calcul_Interpolation<_PARAM_>
285 public : Calcul_Interpolation_Tetra10(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
286 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
288 // ON SUPPOSE IMPLICITEMENT QUE LES NOEUDS SUPPLEMENTAIRES SONT BIEN DES NOEUDS MILIEUX
289 int num0=mailles->DONNE_SOMMET_MAILLE(num_maille,0);
290 int num1=mailles->DONNE_SOMMET_MAILLE(num_maille,1);
291 int num2=mailles->DONNE_SOMMET_MAILLE(num_maille,2);
292 int num3=mailles->DONNE_SOMMET_MAILLE(num_maille,3);
293 int num01=mailles->DONNE_SOMMET_MAILLE(num_maille,4);
294 int num02=mailles->DONNE_SOMMET_MAILLE(num_maille,6);
295 int num03=mailles->DONNE_SOMMET_MAILLE(num_maille,7);
296 int num12=mailles->DONNE_SOMMET_MAILLE(num_maille,5);
297 int num13=mailles->DONNE_SOMMET_MAILLE(num_maille,8);
298 int num23=mailles->DONNE_SOMMET_MAILLE(num_maille,9);
300 double x0=(*noeuds)[num0][0];
301 double y0=(*noeuds)[num0][1];
302 double z0=(*noeuds)[num0][2];
303 double x1=(*noeuds)[num1][0];
304 double y1=(*noeuds)[num1][1];
305 double z1=(*noeuds)[num1][2];
306 double x2=(*noeuds)[num2][0];
307 double y2=(*noeuds)[num2][1];
308 double z2=(*noeuds)[num2][2];
309 double x3=(*noeuds)[num3][0];
310 double y3=(*noeuds)[num3][1];
311 double z3=(*noeuds)[num3][2];
313 VALEURCHAMP v0=(*champ)[num0];
314 VALEURCHAMP v1=(*champ)[num1];
315 VALEURCHAMP v2=(*champ)[num2];
316 VALEURCHAMP v3=(*champ)[num3];
317 VALEURCHAMP v01=(*champ)[num01];
318 VALEURCHAMP v02=(*champ)[num02];
319 VALEURCHAMP v03=(*champ)[num03];
320 VALEURCHAMP v12=(*champ)[num12];
321 VALEURCHAMP v13=(*champ)[num13];
322 VALEURCHAMP v23=(*champ)[num23];
328 double lambda0=(y3*(z1-z2)+y1*(z2-z3)+y2*(z3-z1))*x+(x3*(z2-z1)+x2*(z1-z3)+x1*(z3-z2))*y+(x3*(y1-y2)+x1*(y2-y3)+x2*(y3-y1))*z+((x3*y2-x2*y3)*z1+(x1*y3-x3*y1)*z2+(x2*y1-x1*y2)*z3);
329 double lambda1=-((y0*(z2-z3)+y2*(z3-z0)+y3*(z0-z2))*x+(x0*(z3-z2)+x3*(z2-z0)+x2*(z0-z3))*y+(x0*(y2-y3)+x2*(y3-y0)+x3*(y0-y2))*z+((x0*y3-x2*y0)*z2+(x2*y0-x0*y2)*z3+(x3*y2-x2*y3)*z0));
330 double lambda2=(y1*(z3-z0)+y3*(z0-z1)+y0*(z1-z3))*x+(x1*(z0-z3)+x0*(z3-z1)+x3*(z1-z0))*y+(x1*(y3-y0)+x3*(y0-y1)+x0*(y1-y3))*z+((x1*y0-x2*y1)*z3+(x3*y1-x1*y3)*z0+(x0*y3-x3*y0)*z1);
331 double lambda3=-((y2*(z0-z1)+y0*(z1-z2)+y1*(z2-z0))*x+(x2*(z1-z0)+x1*(z0-z2)+x0*(z2-z1))*y+(x2*(y0-y1)+x0*(y1-y2)+x1*(y2-y0))*z+((x2*y1-x2*y2)*z0+(x0*y2-x2*y0)*z1+(x1*y0-x0*y1)*z2));
334 // A SIMPLFIFIER, FEIGNASSE !
335 double delta = (x1*y2-x1*y3)*z0 +(x0-y3- x0*y2)*z1 + x0*y3*z1 - x1*y0*z2 + x0*y1*z2 - x0*y3*z2 + x1*y3*z2 + x3*(y1*z0 - y2*z0 - y0*z1 + y2*z1 + y0*z2 - y1*z2) + x1*y0*z3 - x0*y1*z3 + x0*y2*z3 - x1*y2*z3 + x2*(y3*z0 + y0*z1 - y3*z1 - y0*z3 + y1*(-z0 + z3));
338 cout<<" ### Pour ( "<<x<<" , "<<y<<" , "<<z<<" )"<<endl;
339 cout<<" ### delta = "<<delta<<endl;
340 cout<<" ### Lambda0 = "<<lambda0<<endl;
341 cout<<" ### Lambda1 = "<<lambda1<<endl;
342 cout<<" ### Lambda2 = "<<lambda2<<endl;
343 cout<<" ### Lambda3 = "<<lambda3<<endl;
346 VALEURCHAMP retour(v0.SIZE()); //
348 retour=2*(lambda0*lambda0*v0+
352 2*(lambda0*lambda1*v01+
357 lambda2*lambda3*v23))/(delta*delta)+
358 (lambda0*v0+lambda1*v1+lambda2*v2+lambda3*v3)/(-delta);
363 _TEMPLATE_ class Calcul_Interpolation_Hexa8 : public Calcul_Interpolation<_PARAM_>
365 public : Calcul_Interpolation_Hexa8(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
366 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
369 int num0=mailles->DONNE_SOMMET_MAILLE(num_maille,0);
370 int num1=mailles->DONNE_SOMMET_MAILLE(num_maille,1);
371 int num2=mailles->DONNE_SOMMET_MAILLE(num_maille,2);
372 int num3=mailles->DONNE_SOMMET_MAILLE(num_maille,3);
373 int num4=mailles->DONNE_SOMMET_MAILLE(num_maille,4);
374 int num5=mailles->DONNE_SOMMET_MAILLE(num_maille,5);
375 int num6=mailles->DONNE_SOMMET_MAILLE(num_maille,6);
376 int num7=mailles->DONNE_SOMMET_MAILLE(num_maille,7);
379 double x0=(*noeuds)[num0][0];double y0=(*noeuds)[num0][1];double z0=(*noeuds)[num0][2];
380 double x1=(*noeuds)[num1][0];double y1=(*noeuds)[num1][1];double z1=(*noeuds)[num1][2];
381 double x2=(*noeuds)[num2][0];double y2=(*noeuds)[num2][1];double z2=(*noeuds)[num2][2];
382 double x3=(*noeuds)[num3][0];double y3=(*noeuds)[num3][1];double z3=(*noeuds)[num3][2];
383 double x4=(*noeuds)[num4][0];double y4=(*noeuds)[num4][1];double z4=(*noeuds)[num4][2];
384 double x5=(*noeuds)[num5][0];double y5=(*noeuds)[num5][1];double z5=(*noeuds)[num5][2];
385 double x6=(*noeuds)[num6][0];double y6=(*noeuds)[num6][1];double z6=(*noeuds)[num6][2];
386 double x7=(*noeuds)[num7][0];double y7=(*noeuds)[num7][1];double z7=(*noeuds)[num7][2];
389 VALEURCHAMP v0=(*champ)[num0];
390 VALEURCHAMP v1=(*champ)[num1];
391 VALEURCHAMP v2=(*champ)[num2];
392 VALEURCHAMP v3=(*champ)[num3];
393 VALEURCHAMP v4=(*champ)[num4];
394 VALEURCHAMP v5=(*champ)[num5];
395 VALEURCHAMP v6=(*champ)[num6];
396 VALEURCHAMP v7=(*champ)[num7];
403 // double face2367(double x,double y,double z) {return (x6*(-y2 + y3) + x3*(y2 - y6) + x2*(-y3 + y6))*z - x6*y3*z2 + x3*y6*z2 + x6*y2*z3 - x2*y6*z3 - x3*y2*z6 + x2*y3*z6 + y*(x6*(z2 - z3) + x2*(z3 - z6) + x3*(-z2 + z6)) + x*(y6*(-z2 + z3) + y3*(z2 - z6) + y2*(-z3 + z6));}
404 // double face4567(double x,double y,double z) {return (x6*(-y4 + y5) + x5*(y4 - y6) + x4*(-y5 + y6))*z - x6*y5*z4 + x5*y6*z4 + x6*y4*z5 - x4*y6*z5 - x5*y4*z6 + x4*y5*z6 + y*(x6*(z4 - z5) + x4*(z5 - z6) + x5*(-z4 + z6)) + x*(y6*(-z4 + z5) + y5*(z4 - z6) + y4*(-z5 + z6));}
405 // double face1256(double x,double y,double z) {return (x5*(-y1 + y2) + x2*(y1 - y5) + x1*(-y2 + y5))*z - x5*y2*z1 + x2*y5*z1 + x5*y1*z2 - x1*y5*z2 - x2*y1*z5 + x1*y2*z5 + y*(x5*(z1 - z2) + x1*(z2 - z5) + x2*(-z1 + z5)) + x*(y5*(-z1 + z2) + y2*(z1 - z5) + y1*(-z2 + z5));}
406 // double face0347(double x,double y,double z) {return (x4*(-y0 + y3) + x3*(y0 - y4) + x0*(-y3 + y4))*z - x4*y3*z0 + x3*y4*z0 + x4*y0*z3 - x0*y4*z3 - x3*y0*z4 + x0*y3*z4 + y*(x4*(z0 - z3) + x0*(z3 - z4) + x3*(-z0 + z4)) + x*(y4*(-z0 + z3) + y3*(z0 - z4) + y0*(-z3 + z4));}
407 // double face0145(double x,double y,double z) {return (x4*(-y0 + y1) + x1*(y0 - y4) + x0*(-y1 + y4))*z - x4*y1*z0 + x1*y4*z0 + x4*y0*z1 - x0*y4*z1 - x1*y0*z4 + x0*y1*z4 + y*(x4*(z0 - z1) + x0*(z1 - z4) + x1*(-z0 + z4)) + x*(y4*(-z0 + z1) + y1*(z0 - z4) + y0*(-z1 + z4));}
408 // double face0123(double x,double y,double z) {return (x2*(-y0 + y1) + x1*(y0 - y2) + x0*(-y1 + y2))*z - x2*y1*z0 + x1*y2*z0 + x2*y0*z1 - x0*y2*z1 - x1*y0*z2 + x0*y1*z2 + y*(x2*(z0 - z1) + x0*(z1 - z2) + x1*(-z0 + z2)) + x*(y2*(-z0 + z1) + y1*(z0 - z2) + y0*(-z1 + z2));}
410 double mu0=face1256(x,y,z)*face2367(x,y,z)*face4567(x,y,z)/(face1256(x0,y0,z0)*face2367(x0,y0,z0)*face4567(x0,y0,z0));
411 double mu1=face0347(x,y,z)*face2367(x,y,z)*face4567(x,y,z)/(face0347(x1,y1,z1)*face2367(x1,y1,z1)*face4567(x1,y1,z1));
412 double mu2=face0347(x,y,z)*face0145(x,y,z)*face4567(x,y,z)/(face0347(x2,y2,z2)*face0145(x2,y2,z2)*face4567(x2,y2,z2));
413 double mu3=face1256(x,y,z)*face0145(x,y,z)*face4567(x,y,z)/(face1256(x3,y3,z3)*face0145(x3,y3,z3)*face4567(x3,y3,z3));
414 double mu4=face1256(x,y,z)*face2367(x,y,z)*face0123(x,y,z)/(face1256(x4,y4,z4)*face2367(x4,y4,z4)*face0123(x4,y4,z4));
415 double mu5=face0347(x,y,z)*face2367(x,y,z)*face0123(x,y,z)/(face0347(x5,y5,z5)*face2367(x5,y5,z5)*face0123(x5,y5,z5));
416 double mu6=face0347(x,y,z)*face0145(x,y,z)*face0123(x,y,z)/(face0347(x6,y6,z6)*face0145(x6,y6,z6)*face0123(x6,y6,z6));
417 double mu7=face1256(x,y,z)*face0145(x,y,z)*face0123(x,y,z)/(face1256(x7,y7,z7)*face0145(x7,y7,z7)*face0123(x7,y7,z7));
419 VALEURCHAMP retour(v0.SIZE()); //
421 retour=(mu0*v0+mu1*v1+mu2*v2+mu3*v3+mu4*v4+mu5*v5+mu6*v6+mu7*v7);
426 _TEMPLATE_ class Calcul_Interpolation_Hexa20 : public Calcul_Interpolation<_PARAM_>
428 public : Calcul_Interpolation_Hexa20(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
429 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
431 cerr<<"Interpolation H2 pasencore implémentée"<<endl;
435 _TEMPLATE_ class Calcul_Interpolation_Penta6 : public Calcul_Interpolation<_PARAM_>
437 public : Calcul_Interpolation_Penta6(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
438 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
440 cerr<<"Interpolation Pe1 pasencore implémentée"<<endl;
444 _TEMPLATE_ class Calcul_Interpolation_Penta15 : public Calcul_Interpolation<_PARAM_>
446 public : Calcul_Interpolation_Penta15(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
447 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
449 cerr<<"Interpolation Pe2 pasencore implémentée"<<endl;
453 _TEMPLATE_ class Calcul_Interpolation_Pyra5 : public Calcul_Interpolation<_PARAM_>
455 public : Calcul_Interpolation_Pyra5(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
456 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
458 cerr<<"Interpolation Py1 pasencore implémentée"<<endl;
462 _TEMPLATE_ class Calcul_Interpolation_Pyra13 : public Calcul_Interpolation<_PARAM_>
464 public : Calcul_Interpolation_Pyra13(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
465 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
467 cerr<<"Interpolation Py2 pasencore implémentée"<<endl;