src/MEDMEM/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 #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)))
  15
  16
  17
  18 // DECLARATION
  19
  20 _TEMPLATE_ class Calcul_Interpolation
  21 {
  22 protected :
  23         NUAGENOEUD * noeuds;
  24         NUAGEMAILLE * mailles;
  25         CHAMP * champ;
  26 public :
  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);}
  31 };
  32
  33 _TEMPLATE_ class Calcul_Interpolation_P0;
  34
  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;
  47
  48 _TEMPLATE_ class Calcul_Hybride
  49 {
  50 protected :
  51         NUAGEMAILLE * mailles;
  52         map<int,Calcul_Interpolation<_PARAM_> *> fonctions;
  53 public :
  54         Calcul_Hybride():mailles(NULL) {}
  55         Calcul_Hybride(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c);
  56         ~Calcul_Hybride() {}
  57         VALEURCHAMP operator() (const NOEUD & n, int num_maille);
  58 };
  59
  60 //CODE
  61
  62 _TEMPLATE_ Calcul_Hybride<_PARAM_>::Calcul_Hybride(NUAGENOEUD * nn,NUAGEMAILLE * nm,CHAMP * c):mailles(nm)
  63         {
  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);
  76         }
  77
  78 _TEMPLATE_ VALEURCHAMP Calcul_Hybride<_PARAM_>::operator() (const NOEUD & n, int num_maille)
  79         {
  80         return fonctions[mailles->DONNE_TYPE_MAILLE(num_maille)]->operator()(n,num_maille);
  81         }
  82
  83 _TEMPLATE_ class Calcul_Interpolation_P0      : public Calcul_Interpolation<_PARAM_>
  84 {
  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)
  87         {
  88         return (*champ)[num_maille];
  89         }
  90 };
  91 _TEMPLATE_ class Calcul_Interpolation_Tria3   : public Calcul_Interpolation<_PARAM_>
  92 {
  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)
  95         {
  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);
  99
 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];
 106
 107         VALEURCHAMP v0=(*champ)[num0];
 108         VALEURCHAMP v1=(*champ)[num1];
 109         VALEURCHAMP v2=(*champ)[num2];
 110
 111         double x=n[0];
 112         double y=n[1];
 113
 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);
 117
 118         double delta = (x2-x1)*y0+(x0-x2)*y1+(x1-x0)*y2;
 119
 120         VALEURCHAMP retour(v0.SIZE());
 121
 122         return (1/delta)*(lambda0*v0+lambda1*v1+lambda2*v2);
 123
 124         // return retour; //
 125         }
 126 };
 127 _TEMPLATE_ class Calcul_Interpolation_Tria6   : public Calcul_Interpolation<_PARAM_>
 128 {
 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)
 131         {
 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);
 139
 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];
 146
 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];
 153
 154         double x=n[0];
 155         double y=n[1];
 156
 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);
 160
 161         double delta = (x2-x1)*y0+(x0-x2)*y1+(x1-x0)*y2;
 162
 163
 164         // VALEURCHAMP retour(v0.SIZE()); //
 165
 166         return  2*(lambda0*lambda0*v0+
 167                    lambda1*lambda1*v1+
 168                    lambda2*lambda2*v2+
 169                    2*(lambda0*lambda1*v01+
 170                       lambda1*lambda2*v12+
 171                       lambda2*lambda0*v20))/(delta*delta)+
 172                 (lambda0*v0+lambda1*v1+lambda2*v2)/(-delta);
 173
 174         // return retour; //
 175         }
 176 };
 177 _TEMPLATE_ class Calcul_Interpolation_Quad4   : public Calcul_Interpolation<_PARAM_>
 178 {
 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)
 181         {
 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);
 186
 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];
 195
 196         VALEURCHAMP v0=(*champ)[num0];
 197         VALEURCHAMP v1=(*champ)[num1];
 198         VALEURCHAMP v2=(*champ)[num2];
 199         VALEURCHAMP v3=(*champ)[num3];
 200
 201         double x=n[0];
 202         double y=n[1];
 203
 204
 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));
 209
 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);
 211
 212         /*
 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;
 219         //*/
 220         // VALEURCHAMP retour(v0.SIZE()); //
 221
 222         return (mu0*v0+mu1*v1+mu2*v2+mu3*v3)/delta;
 223
 224         // return retour; //
 225         }
 226 };
 227 _TEMPLATE_ class Calcul_Interpolation_Quad8   : public Calcul_Interpolation<_PARAM_>
 228 {
 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)
 231         {
 232         cerr<<"Interpolation Q2 pas encore implémentée"<<endl;
 233         exit(-1);
 234         }
 235 };
 236 _TEMPLATE_ class Calcul_Interpolation_Tetra4  : public Calcul_Interpolation<_PARAM_>
 237 {
 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)
 240         {
 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);
 245
 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];
 258
 259         VALEURCHAMP v0=(*champ)[num0];
 260         VALEURCHAMP v1=(*champ)[num1];
 261         VALEURCHAMP v2=(*champ)[num2];
 262         VALEURCHAMP v3=(*champ)[num3];
 263
 264         double x=n[0];
 265         double y=n[1];
 266         double z=n[2];
 267
 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));
 272
 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));
 275
 276         VALEURCHAMP retour(v0.SIZE()); //
 277
 278         retour=(lambda0*v0+lambda1*v1+lambda2*v2+lambda3*v3)/delta;
 279
 280         return retour; //
 281         }
 282 };
 283 _TEMPLATE_ class Calcul_Interpolation_Tetra10 : public Calcul_Interpolation<_PARAM_>
 284 {
 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)
 287         {
 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);
 299
 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];
 312
 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];
 323
 324         double x=n[0];
 325         double y=n[1];
 326         double z=n[2];
 327
 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));
 332
 333
 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));
 336
 337         /*
 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;
 344         //*/
 345
 346         VALEURCHAMP retour(v0.SIZE()); //
 347
 348         retour=2*(lambda0*lambda0*v0+
 349                   lambda1*lambda1*v1+
 350                   lambda2*lambda2*v2+
 351                   lambda3*lambda3*v3+
 352                   2*(lambda0*lambda1*v01+
 353                      lambda0*lambda2*v02+
 354                      lambda0*lambda3*v03+
 355                      lambda1*lambda2*v12+
 356                      lambda1*lambda3*v13+
 357                      lambda2*lambda3*v23))/(delta*delta)+
 358                (lambda0*v0+lambda1*v1+lambda2*v2+lambda3*v3)/(-delta);
 359
 360         return retour; //
 361         }
 362 };
 363 _TEMPLATE_ class Calcul_Interpolation_Hexa8   : public Calcul_Interpolation<_PARAM_>
 364 {
 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)
 367         {
 368
 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);
 377
 378
 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];
 387
 388
 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];
 397
 398         double x=n[0];
 399         double y=n[1];
 400         double z=n[2];
 401
 402
 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));}
 409
 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));
 418
 419         VALEURCHAMP retour(v0.SIZE()); //
 420
 421         retour=(mu0*v0+mu1*v1+mu2*v2+mu3*v3+mu4*v4+mu5*v5+mu6*v6+mu7*v7);
 422
 423         return retour; //
 424         }
 425 };
 426 _TEMPLATE_ class Calcul_Interpolation_Hexa20  : public Calcul_Interpolation<_PARAM_>
 427 {
 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)
 430         {
 431         cerr<<"Interpolation H2 pasencore implémentée"<<endl;
 432         exit(-1);
 433         }
 434 };
 435 _TEMPLATE_ class Calcul_Interpolation_Penta6  : public Calcul_Interpolation<_PARAM_>
 436 {
 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)
 439         {
 440         cerr<<"Interpolation Pe1 pasencore implémentée"<<endl;
 441         exit(-1);
 442         }
 443 };
 444 _TEMPLATE_ class Calcul_Interpolation_Penta15 : public Calcul_Interpolation<_PARAM_>
 445 {
 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)
 448         {
 449         cerr<<"Interpolation Pe2 pasencore implémentée"<<endl;
 450         exit(-1);
 451         }
 452 };
 453 _TEMPLATE_ class Calcul_Interpolation_Pyra5   : public Calcul_Interpolation<_PARAM_>
 454 {
 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)
 457         {
 458         cerr<<"Interpolation Py1 pasencore implémentée"<<endl;
 459         exit(-1);
 460         }
 461 };
 462 _TEMPLATE_ class Calcul_Interpolation_Pyra13  : public Calcul_Interpolation<_PARAM_>
 463 {
 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)
 466         {
 467         cerr<<"Interpolation Py2 pasencore implémentée"<<endl;
 468         exit(-1);
 469         }
 470 };
 471
 472 #undef _TEMPLATE_
 473 #undef _PARAM_
 474
 475 #endif