src/INTERPOLATION/MEDMEM_InterpolationTools.hxx

   1 #ifndef MEDMEM_INTERPOLATION_TOOLS_HXX
   2 #define MEDMEM_INTERPOLATION_TOOLS_HXX
   3
   4 #include "MEDMEM_define.hxx"
   5
   6 #define _TEMPLATE_ template <class CHAMP,class VALEURCHAMP,class NUAGENOEUD, class NOEUD, class NUAGEMAILLE>
   7 #define _PARAM_ CHAMP,VALEURCHAMP,NUAGENOEUD,NOEUD,NUAGEMAILLE
   8
   9 // ces macros définissent pour une face carrée plane la fonction projection sur cette face, non normalisée, la numérotation est dans le cas d'un hexaedre
  10 #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)))
  11 #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)))
  12 #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)))
  13 #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)))
  14 #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)))
  15 #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)))
  16 // des macros définissent pour une face triangulaire orientée vers l'extérieur de la maille la fonction de projection, non normalisée ( =(12^13).1M )
  17 #define face(x1,y1,z1,x2,y2,z2,x3,y3,z3,x,y,z)  ( ((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1))*(x-x1)+((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1))*(y-y1)+((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1))*(z-z1) )
  18 #define projection(x1,y1,z1,x2,y2,z2,x3,y3,z3,x0,y0,z0,x,y,z) (face(x1,y1,z1,x2,y2,z2,x3,y3,z3,x,y,z)/face(x1,y1,z1,x2,y2,z2,x3,y3,z3,x0,y0,z0))
  19
  20 // DECLARATION
  21
  22 _TEMPLATE_ class Calcul_Interpolation
  23 {
  24 protected :
  25         NUAGENOEUD * noeuds;
  26         NUAGEMAILLE * mailles;
  27         CHAMP * champ;
  28 public :
  29         Calcul_Interpolation():noeuds(NULL),mailles(NULL),champ(NULL) {}
  30         Calcul_Interpolation(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):noeuds(nn),mailles(nm),champ(c) {}
  31         ~Calcul_Interpolation() {}
  32         virtual VALEURCHAMP operator() (const NOEUD & n, int num_maille) {cerr<<"APPEL OPERATOR() DE LA CLASSE MERE CALCUL_INTERPOLATION => EXIT(-1)"<<endl;exit(-1);}
  33 };
  34
  35 _TEMPLATE_ class Calcul_Interpolation_P0;
  36
  37 _TEMPLATE_ class Calcul_Interpolation_Tria3;
  38 _TEMPLATE_ class Calcul_Interpolation_Tria6;
  39 _TEMPLATE_ class Calcul_Interpolation_Quad4;
  40 _TEMPLATE_ class Calcul_Interpolation_Quad8;
  41 _TEMPLATE_ class Calcul_Interpolation_Tetra4;
  42 _TEMPLATE_ class Calcul_Interpolation_Tetra10;
  43 _TEMPLATE_ class Calcul_Interpolation_Hexa8;
  44 _TEMPLATE_ class Calcul_Interpolation_Hexa20;
  45 _TEMPLATE_ class Calcul_Interpolation_Penta6;
  46 _TEMPLATE_ class Calcul_Interpolation_Penta15;
  47 _TEMPLATE_ class Calcul_Interpolation_Pyra5;
  48 _TEMPLATE_ class Calcul_Interpolation_Pyra13;
  49
  50 _TEMPLATE_ class Calcul_Hybride
  51 {
  52 protected :
  53         NUAGEMAILLE * mailles;
  54         map<int,Calcul_Interpolation<_PARAM_> *> fonctions;
  55 public :
  56         Calcul_Hybride():mailles(NULL) {}
  57         Calcul_Hybride(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c);
  58         ~Calcul_Hybride() {}
  59         VALEURCHAMP operator() (const NOEUD & n, int num_maille);
  60 };
  61
  62 //CODE
  63
  64 _TEMPLATE_ Calcul_Hybride<_PARAM_>::Calcul_Hybride(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):mailles(nm)
  65         {
  66         fonctions[ MED_EN::MED_TRIA3  ]=new Calcul_Interpolation_Tria3  <_PARAM_>(nn,nm,c);
  67         fonctions[ MED_EN::MED_TRIA6  ]=new Calcul_Interpolation_Tria6  <_PARAM_>(nn,nm,c);
  68         fonctions[ MED_EN::MED_QUAD4  ]=new Calcul_Interpolation_Quad4  <_PARAM_>(nn,nm,c);
  69         fonctions[ MED_EN::MED_QUAD8  ]=new Calcul_Interpolation_Quad8  <_PARAM_>(nn,nm,c);
  70         fonctions[ MED_EN::MED_TETRA4 ]=new Calcul_Interpolation_Tetra4 <_PARAM_>(nn,nm,c);
  71         fonctions[ MED_EN::MED_TETRA10]=new Calcul_Interpolation_Tetra10<_PARAM_>(nn,nm,c);
  72         fonctions[ MED_EN::MED_HEXA8  ]=new Calcul_Interpolation_Hexa8  <_PARAM_>(nn,nm,c);
  73         fonctions[ MED_EN::MED_HEXA20 ]=new Calcul_Interpolation_Hexa20 <_PARAM_>(nn,nm,c);
  74         fonctions[ MED_EN::MED_PENTA6 ]=new Calcul_Interpolation_Penta6 <_PARAM_>(nn,nm,c);
  75         fonctions[ MED_EN::MED_PENTA15]=new Calcul_Interpolation_Penta15<_PARAM_>(nn,nm,c);
  76         fonctions[ MED_EN::MED_PYRA5  ]=new Calcul_Interpolation_Pyra5  <_PARAM_>(nn,nm,c);
  77         fonctions[ MED_EN::MED_PYRA13 ]=new Calcul_Interpolation_Pyra13 <_PARAM_>(nn,nm,c);
  78         }
  79
  80 _TEMPLATE_ VALEURCHAMP Calcul_Hybride<_PARAM_>::operator() (const NOEUD & n, int num_maille)
  81         {
  82         return fonctions[(*mailles)[num_maille].DONNE_TYPE_MED_MAILLE()]->operator()(n,num_maille);
  83         }
  84
  85 _TEMPLATE_ class Calcul_Interpolation_P0      : public Calcul_Interpolation<_PARAM_>
  86 {
  87 public : Calcul_Interpolation_P0(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
  88 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
  89         {
  90         return (*Calcul_Interpolation<_PARAM_>::champ)[num_maille];
  91         }
  92 };
  93 _TEMPLATE_ class Calcul_Interpolation_Tria3   : public Calcul_Interpolation<_PARAM_>
  94 {
  95 public : Calcul_Interpolation_Tria3(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
  96 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
  97         {
  98         int num0=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][0];
  99         int num1=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][1];
 100         int num2=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][2];
 101
 102         double x0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][0];
 103         double y0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][1];
 104         double x1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][0];
 105         double y1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][1];
 106         double x2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][0];
 107         double y2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][1];
 108
 109         VALEURCHAMP v0=(*Calcul_Interpolation<_PARAM_>::champ)[num0];
 110         VALEURCHAMP v1=(*Calcul_Interpolation<_PARAM_>::champ)[num1];
 111         VALEURCHAMP v2=(*Calcul_Interpolation<_PARAM_>::champ)[num2];
 112
 113         double x=n[0];
 114         double y=n[1];
 115
 116         double lambda0=(y1-y2)*x+(x2-x1)*y+(x1*y2-x2*y1);
 117         double lambda1=(y2-y0)*x+(x0-x2)*y+(x2*y0-x0*y2);
 118         double lambda2=(y0-y1)*x+(x1-x0)*y+(x0*y1-x1*y0);
 119
 120         double delta = (x2-x1)*y0+(x0-x2)*y1+(x1-x0)*y2;
 121
 122         VALEURCHAMP retour(v0.SIZE());
 123
 124         retour=(1/delta)*(lambda0*v0+lambda1*v1+lambda2*v2);
 125
 126         return retour; //
 127         }
 128 };
 129 _TEMPLATE_ class Calcul_Interpolation_Tria6   : public Calcul_Interpolation<_PARAM_>
 130 {
 131 public : Calcul_Interpolation_Tria6(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
 132 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
 133         {
 134         // ON SUPPOSE IMPLICITEMENT QUE LES NOEUDS SUPPLEMENTAIRES SONT BIEN DES NOEUDS MILIEUX
 135         int num0 =(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][0];
 136         int num1 =(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][1];
 137         int num2 =(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][2];
 138         int num01=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][3];
 139         int num12=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][4];
 140         int num20=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][5];
 141
 142         double x0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][0];
 143         double y0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][1];
 144         double x1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][0];
 145         double y1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][1];
 146         double x2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][0];
 147         double y2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][1];
 148
 149         VALEURCHAMP v0=(*Calcul_Interpolation<_PARAM_>::champ)[num0];
 150         VALEURCHAMP v1=(*Calcul_Interpolation<_PARAM_>::champ)[num1];
 151         VALEURCHAMP v2=(*Calcul_Interpolation<_PARAM_>::champ)[num2];
 152         VALEURCHAMP v01=(*Calcul_Interpolation<_PARAM_>::champ)[num01];
 153         VALEURCHAMP v12=(*Calcul_Interpolation<_PARAM_>::champ)[num12];
 154         VALEURCHAMP v20=(*Calcul_Interpolation<_PARAM_>::champ)[num20];
 155
 156         double x=n[0];
 157         double y=n[1];
 158
 159         double lambda0=(y1-y2)*x+(x2-x1)*y+(x1*y2-x2*y1);
 160         double lambda1=(y2-y0)*x+(x0-x2)*y+(x2*y0-x0*y2);
 161         double lambda2=(y0-y1)*x+(x1-x0)*y+(x0*y1-x1*y0);
 162
 163         double delta = (x2-x1)*y0+(x0-x2)*y1+(x1-x0)*y2;
 164
 165
 166         // VALEURCHAMP retour(v0.SIZE()); //
 167
 168         return  2*(lambda0*lambda0*v0+
 169                    lambda1*lambda1*v1+
 170                    lambda2*lambda2*v2+
 171                    2*(lambda0*lambda1*v01+
 172                       lambda1*lambda2*v12+
 173                       lambda2*lambda0*v20))/(delta*delta)+
 174                 (lambda0*v0+lambda1*v1+lambda2*v2)/(-delta);
 175
 176         // return retour; //
 177         }
 178 };
 179 _TEMPLATE_ class Calcul_Interpolation_Quad4   : public Calcul_Interpolation<_PARAM_>
 180 {
 181 public : Calcul_Interpolation_Quad4(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
 182 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
 183         {
 184         int num0=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][0];
 185         int num1=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][1];
 186         int num2=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][2];
 187         int num3=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][3];
 188
 189         double x0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][0];
 190         double y0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][1];
 191         double x1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][0];
 192         double y1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][1];
 193         double x2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][0];
 194         double y2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][1];
 195         double x3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][0];
 196         double y3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][1];
 197
 198         VALEURCHAMP v0=(*Calcul_Interpolation<_PARAM_>::champ)[num0];
 199         VALEURCHAMP v1=(*Calcul_Interpolation<_PARAM_>::champ)[num1];
 200         VALEURCHAMP v2=(*Calcul_Interpolation<_PARAM_>::champ)[num2];
 201         VALEURCHAMP v3=(*Calcul_Interpolation<_PARAM_>::champ)[num3];
 202
 203         double x=n[0];
 204         double y=n[1];
 205
 206
 207         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)));
 208         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));
 209         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)));
 210         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));
 211
 212         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
 214         /*
 215         cout<<"  ### Pour ( "<<x<<" , "<<y<<" )"<<endl;
 216         cout<<"  ### delta = "<<delta<<endl;
 217         cout<<"  ### Mu0 = "<<mu0<<endl;
 218         cout<<"  ### Mu1 = "<<mu1<<endl;
 219         cout<<"  ### Mu2 = "<<mu2<<endl;
 220         cout<<"  ### Mu3 = "<<mu3<<endl;
 221         */
 222         VALEURCHAMP retour(v0.SIZE()); //
 223
 224         retour=(mu0*v0+mu1*v1+mu2*v2+mu3*v3)/delta;
 225
 226         return retour; //
 227         }
 228 };
 229 _TEMPLATE_ class Calcul_Interpolation_Quad8   : public Calcul_Interpolation<_PARAM_>
 230 {
 231 public : Calcul_Interpolation_Quad8(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
 232 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
 233         {
 234         cerr<<"Interpolation Q2 pas encore implémentée"<<endl;
 235         exit(-1);
 236         }
 237 };
 238 _TEMPLATE_ class Calcul_Interpolation_Tetra4  : public Calcul_Interpolation<_PARAM_>
 239 {
 240 public : Calcul_Interpolation_Tetra4(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
 241 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
 242         {
 243         int num0=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][0];
 244         int num1=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][1];
 245         int num2=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][2];
 246         int num3=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][3];
 247
 248         double x0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][0];
 249         double y0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][1];
 250         double z0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][2];
 251         double x1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][0];
 252         double y1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][1];
 253         double z1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][2];
 254         double x2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][0];
 255         double y2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][1];
 256         double z2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][2];
 257         double x3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][0];
 258         double y3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][1];
 259         double z3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][2];
 260
 261         VALEURCHAMP v0=(*Calcul_Interpolation<_PARAM_>::champ)[num0];
 262         VALEURCHAMP v1=(*Calcul_Interpolation<_PARAM_>::champ)[num1];
 263         VALEURCHAMP v2=(*Calcul_Interpolation<_PARAM_>::champ)[num2];
 264         VALEURCHAMP v3=(*Calcul_Interpolation<_PARAM_>::champ)[num3];
 265
 266         double x=n[0];
 267         double y=n[1];
 268         double z=n[2];
 269
 270         double lambda0=face(x1,y1,z1,x2,y2,z2,x3,y3,z3,x,y,z)/face(x1,y1,z1,x2,y2,z2,x3,y3,z3,x0,y0,z0);
 271         double lambda1=face(x0,y0,z0,x2,y2,z2,x3,y3,z3,x,y,z)/face(x0,y0,z0,x2,y2,z2,x3,y3,z3,x1,y1,z1);
 272         double lambda2=face(x1,y1,z1,x0,y0,z0,x3,y3,z3,x,y,z)/face(x1,y1,z1,x0,y0,z0,x3,y3,z3,x2,y2,z2);
 273         double lambda3=face(x1,y1,z1,x2,y2,z2,x0,y0,z0,x,y,z)/face(x1,y1,z1,x2,y2,z2,x0,y0,z0,x3,y3,z3);
 274
 275         VALEURCHAMP retour(v0.SIZE()); //
 276
 277         retour=(lambda0*v0+lambda1*v1+lambda2*v2+lambda3*v3);
 278
 279         return retour; //
 280         }
 281 };
 282 _TEMPLATE_ class Calcul_Interpolation_Tetra10 : public Calcul_Interpolation<_PARAM_>
 283 {
 284 public : Calcul_Interpolation_Tetra10(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
 285 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
 286         {
 287         // ON SUPPOSE IMPLICITEMENT QUE LES NOEUDS SUPPLEMENTAIRES SONT BIEN DES NOEUDS MILIEUX
 288         int num0 =(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][0];
 289         int num1 =(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][1];
 290         int num2 =(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][2];
 291         int num3 =(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][3];
 292         int num01=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][4];
 293         int num02=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][6];
 294         int num03=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][7];
 295         int num12=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][5];
 296         int num13=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][8];
 297         int num23=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][9];
 298
 299         double x0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][0];double y0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][1];double z0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][2];
 300         double x1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][0];double y1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][1];double z1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][2];
 301         double x2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][0];double y2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][1];double z2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][2];
 302         double x3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][0];double y3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][1];double z3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][2];
 303
 304         VALEURCHAMP v0=(*Calcul_Interpolation<_PARAM_>::champ)[num0];
 305         VALEURCHAMP v1=(*Calcul_Interpolation<_PARAM_>::champ)[num1];
 306         VALEURCHAMP v2=(*Calcul_Interpolation<_PARAM_>::champ)[num2];
 307         VALEURCHAMP v3=(*Calcul_Interpolation<_PARAM_>::champ)[num3];
 308         VALEURCHAMP v01=(*Calcul_Interpolation<_PARAM_>::champ)[num01];
 309         VALEURCHAMP v02=(*Calcul_Interpolation<_PARAM_>::champ)[num02];
 310         VALEURCHAMP v03=(*Calcul_Interpolation<_PARAM_>::champ)[num03];
 311         VALEURCHAMP v12=(*Calcul_Interpolation<_PARAM_>::champ)[num12];
 312         VALEURCHAMP v13=(*Calcul_Interpolation<_PARAM_>::champ)[num13];
 313         VALEURCHAMP v23=(*Calcul_Interpolation<_PARAM_>::champ)[num23];
 314
 315         double x=n[0];
 316         double y=n[1];
 317         double z=n[2];
 318
 319         double lambda0=face(x1,y1,z1,x2,y2,z2,x3,y3,z3,x,y,z)/face(x1,y1,z1,x2,y2,z2,x3,y3,z3,x0,y0,z0);
 320         double lambda1=face(x0,y0,z0,x2,y2,z2,x3,y3,z3,x,y,z)/face(x0,y0,z0,x2,y2,z2,x3,y3,z3,x1,y1,z1);
 321         double lambda2=face(x1,y1,z1,x0,y0,z0,x3,y3,z3,x,y,z)/face(x1,y1,z1,x0,y0,z0,x3,y3,z3,x2,y2,z2);
 322         double lambda3=face(x1,y1,z1,x2,y2,z2,x0,y0,z0,x,y,z)/face(x1,y1,z1,x2,y2,z2,x0,y0,z0,x3,y3,z3);
 323
 324         /*
 325         cout<<"  ### Pour ( "<<x<<" , "<<y<<" , "<<z<<" )"<<endl;
 326         cout<<"  ### Lambda0 = "<<lambda0<<endl;
 327         cout<<"  ### Lambda1 = "<<lambda1<<endl;
 328         cout<<"  ### Lambda2 = "<<lambda2<<endl;
 329         cout<<"  ### Lambda3 = "<<lambda3<<endl;
 330         cout<<"  ### Controle = "<<(lambda0+lambda1+lambda2+lambda3)<<endl;
 331         */
 332
 333         VALEURCHAMP retour(v0.SIZE()); //
 334
 335         retour=2*(lambda0*lambda0*v0+
 336                   lambda1*lambda1*v1+
 337                   lambda2*lambda2*v2+
 338                   lambda3*lambda3*v3+
 339                   2*(lambda0*lambda1*v01+
 340                      lambda0*lambda2*v02+
 341                      lambda0*lambda3*v03+
 342                      lambda1*lambda2*v12+
 343                      lambda1*lambda3*v13+
 344                      lambda2*lambda3*v23))+
 345                (-1.0)*(lambda0*v0+lambda1*v1+lambda2*v2+lambda3*v3);
 346
 347         return retour; //
 348         }
 349 };
 350 _TEMPLATE_ class Calcul_Interpolation_Hexa8   : public Calcul_Interpolation<_PARAM_>
 351 {
 352 public : Calcul_Interpolation_Hexa8(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
 353 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
 354         {
 355
 356         // Tres probablement numériquement mauvais, à revoir
 357
 358         int num0=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][0];
 359         int num1=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][1];
 360         int num2=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][2];
 361         int num3=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][3];
 362         int num4=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][4];
 363         int num5=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][5];
 364         int num6=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][6];
 365         int num7=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][7];
 366
 367
 368         double x0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][0];double y0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][1];double z0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][2];
 369         double x1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][0];double y1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][1];double z1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][2];
 370         double x2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][0];double y2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][1];double z2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][2];
 371         double x3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][0];double y3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][1];double z3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][2];
 372         double x4=(*Calcul_Interpolation<_PARAM_>::noeuds)[num4][0];double y4=(*Calcul_Interpolation<_PARAM_>::noeuds)[num4][1];double z4=(*Calcul_Interpolation<_PARAM_>::noeuds)[num4][2];
 373         double x5=(*Calcul_Interpolation<_PARAM_>::noeuds)[num5][0];double y5=(*Calcul_Interpolation<_PARAM_>::noeuds)[num5][1];double z5=(*Calcul_Interpolation<_PARAM_>::noeuds)[num5][2];
 374         double x6=(*Calcul_Interpolation<_PARAM_>::noeuds)[num6][0];double y6=(*Calcul_Interpolation<_PARAM_>::noeuds)[num6][1];double z6=(*Calcul_Interpolation<_PARAM_>::noeuds)[num6][2];
 375         double x7=(*Calcul_Interpolation<_PARAM_>::noeuds)[num7][0];double y7=(*Calcul_Interpolation<_PARAM_>::noeuds)[num7][1];double z7=(*Calcul_Interpolation<_PARAM_>::noeuds)[num7][2];
 376
 377
 378         VALEURCHAMP v0=(*Calcul_Interpolation<_PARAM_>::champ)[num0];
 379         VALEURCHAMP v1=(*Calcul_Interpolation<_PARAM_>::champ)[num1];
 380         VALEURCHAMP v2=(*Calcul_Interpolation<_PARAM_>::champ)[num2];
 381         VALEURCHAMP v3=(*Calcul_Interpolation<_PARAM_>::champ)[num3];
 382         VALEURCHAMP v4=(*Calcul_Interpolation<_PARAM_>::champ)[num4];
 383         VALEURCHAMP v5=(*Calcul_Interpolation<_PARAM_>::champ)[num5];
 384         VALEURCHAMP v6=(*Calcul_Interpolation<_PARAM_>::champ)[num6];
 385         VALEURCHAMP v7=(*Calcul_Interpolation<_PARAM_>::champ)[num7];
 386
 387         double x=n[0];
 388         double y=n[1];
 389         double z=n[2];
 390
 391         // 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));
 392         // 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));
 393         // 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));
 394         // 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));
 395         // 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));
 396         // 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));
 397         // 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));
 398         // 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));
 399
 400         double mu0=projection(x1,y1,z1,x2,y2,z2,x6,y6,z6,x0,y0,z0,x,y,z)*projection(x2,y2,z2,x3,y3,z3,x7,y7,z7,x0,y0,z0,x,y,z)*projection(x4,y4,z4,x5,y5,z5,x6,y6,z6,x0,y0,z0,x,y,z);
 401         double mu1=projection(x0,y0,z0,x3,y3,z3,x7,y7,z7,x1,y1,z1,x,y,z)*projection(x2,y2,z2,x3,y3,z3,x7,y7,z7,x1,y1,z1,x,y,z)*projection(x4,y4,z4,x5,y5,z5,x6,y6,z6,x1,y1,z1,x,y,z);
 402         double mu2=projection(x0,y0,z0,x3,y3,z3,x7,y7,z7,x2,y2,z2,x,y,z)*projection(x0,y0,z0,x1,y1,z1,x5,y5,z5,x2,y2,z2,x,y,z)*projection(x4,y4,z4,x5,y5,z5,x6,y6,z6,x2,y2,z2,x,y,z);
 403         double mu3=projection(x1,y1,z1,x2,y2,z2,x6,y6,z6,x3,y3,z3,x,y,z)*projection(x0,y0,z0,x1,y1,z1,x5,y5,z5,x3,y3,z3,x,y,z)*projection(x4,y4,z4,x6,y6,z6,x7,y7,z7,x3,y3,z3,x,y,z);
 404         double mu4=projection(x2,y2,z2,x6,y6,z6,x5,y5,z5,x4,y4,z4,x,y,z)*projection(x2,y2,z2,x3,y3,z3,x7,y7,z7,x4,y4,z4,x,y,z)*projection(x0,y0,z0,x1,y1,z1,x2,y2,z2,x4,y4,z4,x,y,z);
 405         double mu5=projection(x3,y3,z3,x7,y7,z7,x4,y4,z4,x5,y5,z5,x,y,z)*projection(x2,y2,z2,x3,y3,z3,x6,y6,z6,x5,y5,z5,x,y,z)*projection(x0,y0,z0,x1,y1,z1,x2,y2,z2,x5,y5,z5,x,y,z);
 406         double mu6=projection(x3,y3,z3,x7,y7,z7,x4,y4,z4,x6,y6,z6,x,y,z)*projection(x1,y1,z1,x5,y5,z5,x4,y4,z4,x6,y6,z6,x,y,z)*projection(x0,y0,z0,x1,y1,z1,x2,y2,z2,x6,y6,z6,x,y,z);
 407         double mu7=projection(x2,y2,z2,x6,y6,z6,x5,y5,z5,x7,y7,z7,x,y,z)*projection(x1,y1,z1,x5,y5,z5,x4,y4,z4,x7,y7,z7,x,y,z)*projection(x0,y0,z0,x1,y1,z1,x3,y3,z3,x7,y7,z7,x,y,z);
 408
 409         VALEURCHAMP retour(v0.SIZE()); //
 410
 411         retour=(mu0*v0+mu1*v1+mu2*v2+mu3*v3+mu4*v4+mu5*v5+mu6*v6+mu7*v7);
 412
 413         return retour; //
 414         }
 415 };
 416 _TEMPLATE_ class Calcul_Interpolation_Hexa20  : public Calcul_Interpolation<_PARAM_>
 417 {
 418 public : Calcul_Interpolation_Hexa20(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
 419 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
 420         {
 421         cerr<<"Interpolation H2 pasencore implémentée"<<endl;
 422         exit(-1);
 423         }
 424 };
 425 _TEMPLATE_ class Calcul_Interpolation_Penta6  : public Calcul_Interpolation<_PARAM_>
 426 {
 427 public : Calcul_Interpolation_Penta6(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
 428 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
 429         {
 430         int num0=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][0];
 431         int num1=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][1];
 432         int num2=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][2];
 433         int num3=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][3];
 434         int num4=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][4];
 435         int num5=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][5];
 436
 437         double x0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][0];double y0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][1];double z0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][2];
 438         double x1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][0];double y1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][1];double z1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][2];
 439         double x2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][0];double y2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][1];double z2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][2];
 440         double x3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][0];double y3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][1];double z3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][2];
 441         double x4=(*Calcul_Interpolation<_PARAM_>::noeuds)[num4][0];double y4=(*Calcul_Interpolation<_PARAM_>::noeuds)[num4][1];double z4=(*Calcul_Interpolation<_PARAM_>::noeuds)[num4][2];
 442         double x5=(*Calcul_Interpolation<_PARAM_>::noeuds)[num5][0];double y5=(*Calcul_Interpolation<_PARAM_>::noeuds)[num5][1];double z5=(*Calcul_Interpolation<_PARAM_>::noeuds)[num5][2];
 443
 444         VALEURCHAMP v0=(*Calcul_Interpolation<_PARAM_>::champ)[num0];
 445         VALEURCHAMP v1=(*Calcul_Interpolation<_PARAM_>::champ)[num1];
 446         VALEURCHAMP v2=(*Calcul_Interpolation<_PARAM_>::champ)[num2];
 447         VALEURCHAMP v3=(*Calcul_Interpolation<_PARAM_>::champ)[num3];
 448         VALEURCHAMP v4=(*Calcul_Interpolation<_PARAM_>::champ)[num4];
 449         VALEURCHAMP v5=(*Calcul_Interpolation<_PARAM_>::champ)[num5];
 450
 451         double x=n[0];
 452         double y=n[1];
 453         double z=n[2];
 454
 455         double mu0=face(x1,y1,z1,x2,y2,z2,x5,y5,z5,x,y,z)*face(x3,y3,z3,x4,y4,z4,x5,y5,z5,x,y,z)/(face(x1,y1,z1,x2,y2,z2,x5,y5,z5,x0,y0,z0)*face(x3,y3,z3,x4,y4,z4,x5,y5,z5,x0,y0,z0));
 456         double mu1=face(x0,y0,z0,x2,y2,z2,x5,y5,z5,x,y,z)*face(x3,y3,z3,x4,y4,z4,x5,y5,z5,x,y,z)/(face(x0,y0,z0,x2,y2,z2,x5,y5,z5,x1,y1,z1)*face(x3,y3,z3,x4,y4,z4,x5,y5,z5,x1,y1,z1));
 457         double mu2=face(x0,y0,z0,x1,y1,z1,x4,y4,z4,x,y,z)*face(x3,y3,z3,x4,y4,z4,x5,y5,z5,x,y,z)/(face(x0,y0,z0,x1,y1,z1,x4,y4,z4,x2,y2,z2)*face(x3,y3,z3,x4,y4,z4,x5,y5,z5,x2,y2,z2));
 458         double mu3=face(x1,y1,z1,x2,y2,z2,x5,y5,z5,x,y,z)*face(x0,y0,z0,x1,y1,z1,x2,y2,z2,x,y,z)/(face(x1,y1,z1,x2,y2,z2,x5,y5,z5,x3,y3,z3)*face(x0,y0,z0,x1,y1,z1,x2,y2,z2,x3,y3,z3));
 459         double mu4=face(x0,y0,z0,x2,y2,z2,x5,y5,z5,x,y,z)*face(x0,y0,z0,x1,y1,z1,x2,y2,z2,x,y,z)/(face(x0,y0,z0,x2,y2,z2,x5,y5,z5,x4,y4,z4)*face(x0,y0,z0,x1,y1,z1,x2,y2,z2,x4,y4,z4));
 460         double mu5=face(x0,y0,z0,x1,y1,z1,x4,y4,z4,x,y,z)*face(x0,y0,z0,x1,y1,z1,x2,y2,z2,x,y,z)/(face(x0,y0,z0,x1,y1,z1,x4,y4,z4,x5,y5,z5)*face(x0,y0,z0,x1,y1,z1,x2,y2,z2,x5,y5,z5));
 461
 462         VALEURCHAMP retour(v0.SIZE()); //
 463
 464         retour=(mu0*v0+mu1*v1+mu2*v2+mu3*v3+mu4*v4+mu5*v5);
 465
 466         return retour; //
 467         }
 468 };
 469 _TEMPLATE_ class Calcul_Interpolation_Penta15 : public Calcul_Interpolation<_PARAM_>
 470 {
 471 public : Calcul_Interpolation_Penta15(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
 472 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
 473         {
 474         cerr<<"Interpolation Pe2 pasencore implémentée"<<endl;
 475         exit(-1);
 476         }
 477 };
 478 _TEMPLATE_ class Calcul_Interpolation_Pyra5   : public Calcul_Interpolation<_PARAM_>
 479 {
 480 public : Calcul_Interpolation_Pyra5(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
 481 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
 482         {
 483         // NON TESTE
 484         int num0=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][0];
 485         int num1=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][1];
 486         int num2=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][2];
 487         int num3=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][3];
 488         int num4=(*Calcul_Interpolation<_PARAM_>::mailles)[num_maille][4];
 489
 490         double x0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][0];double y0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][1];double z0=(*Calcul_Interpolation<_PARAM_>::noeuds)[num0][2];
 491         double x1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][0];double y1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][1];double z1=(*Calcul_Interpolation<_PARAM_>::noeuds)[num1][2];
 492         double x2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][0];double y2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][1];double z2=(*Calcul_Interpolation<_PARAM_>::noeuds)[num2][2];
 493         double x3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][0];double y3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][1];double z3=(*Calcul_Interpolation<_PARAM_>::noeuds)[num3][2];
 494         double x4=(*Calcul_Interpolation<_PARAM_>::noeuds)[num4][0];double y4=(*Calcul_Interpolation<_PARAM_>::noeuds)[num4][1];double z4=(*Calcul_Interpolation<_PARAM_>::noeuds)[num4][2];
 495
 496         VALEURCHAMP v0=(*Calcul_Interpolation<_PARAM_>::champ)[num0];
 497         VALEURCHAMP v1=(*Calcul_Interpolation<_PARAM_>::champ)[num1];
 498         VALEURCHAMP v2=(*Calcul_Interpolation<_PARAM_>::champ)[num2];
 499         VALEURCHAMP v3=(*Calcul_Interpolation<_PARAM_>::champ)[num3];
 500         VALEURCHAMP v4=(*Calcul_Interpolation<_PARAM_>::champ)[num4];
 501
 502         double x=n[0];
 503         double y=n[1];
 504         double z=n[2];
 505
 506         double mu0=face(x1,y1,z1,x2,y2,z2,x4,y4,z4,x,y,z)*face(x2,y2,z2,x3,y3,z3,x4,y4,z4,x,y,z)/(face(x1,y1,z1,x2,y2,z2,x4,y4,z4,x0,y0,z0)*face(x2,y2,z2,x3,y3,z3,x4,y4,z4,x0,y0,z0));
 507         double mu1=face(x0,y0,z0,x3,y3,z3,x4,y4,z4,x,y,z)*face(x2,y2,z2,x3,y3,z3,x4,y4,z4,x,y,z)/(face(x0,y0,z0,x3,y3,z3,x4,y4,z4,x1,y1,z1)*face(x2,y2,z2,x3,y3,z3,x4,y4,z4,x1,y1,z1));
 508         double mu2=face(x0,y0,z0,x1,y1,z1,x4,y4,z4,x,y,z)*face(x0,y0,z0,x3,y3,z3,x4,y4,z4,x,y,z)/(face(x0,y0,z0,x1,y1,z1,x4,y4,z4,x2,y2,z2)*face(x0,y0,z0,x3,y3,z3,x4,y4,z4,x2,y2,z2));
 509         double mu3=face(x0,y0,z0,x1,y1,z1,x4,y4,z4,x,y,z)*face(x1,y1,z1,x2,y2,z2,x4,y4,z4,x,y,z)/(face(x0,y0,z0,x1,y1,z1,x4,y4,z4,x3,y3,z3)*face(x1,y1,z1,x2,y2,z2,x4,y4,z4,x3,y3,z3));
 510         double mu4=face(x0,y0,z0,x1,y1,z1,x2,y2,z2,x,y,z)/face(x0,y0,z0,x1,y1,z1,x2,y2,z2,x4,y4,z4);
 511
 512         VALEURCHAMP retour(v0.SIZE()); //
 513
 514         retour=(mu0*v0+mu1*v1+mu2*v2+mu3*v3+mu4*v4);
 515
 516         return retour; //
 517         }
 518 };
 519 _TEMPLATE_ class Calcul_Interpolation_Pyra13  : public Calcul_Interpolation<_PARAM_>
 520 {
 521 public : Calcul_Interpolation_Pyra13(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):Calcul_Interpolation<_PARAM_>(nn,nm,c) {}
 522 public : VALEURCHAMP operator() (const NOEUD & n, int num_maille)
 523         {
 524         cerr<<"Interpolation Py2 pasencore implémentée"<<endl;
 525         exit(-1);
 526         }
 527 };
 528
 529 #undef _TEMPLATE_
 530 #undef _PARAM_
 531
 532 #endif