src/MEDMEM/INTERPOLATION/MEDMEM_InterpolationHighLevelObjects.hxx

   1 #ifndef MEDMEM_INTERPOLATION_HIGHLEVEL_OBJECTS_HXX
   2 #define MEDMEM_INTERPOLATION_HIGHLEVEL_OBJECTS_HXX
   3
   4 #include "MEDMEM_Connectivity.hxx"
   5 #include "MEDMEM_WrapperConnectivity.hxx"
   6 #include "MEDMEM_dTree.hxx"
   7 #include "MEDMEM_WrapperNodes.hxx"
   8 #include "MEDMEM_WrapperMesh.hxx"
   9 #include "MEDMEM_WrapperCells.hxx"
  10 #include "MEDMEM_Mapping.hxx"
  11 #include "MEDMEM_WrapperField.hxx"
  12 #include "MEDMEM_InterpolationTools.hxx"
  13
  14 //////////////////////////////////////////////////////////////////
  15 ///                                                            ///
  16 ///                        DECLARATIONS                        ///
  17 ///                                                            ///
  18 //////////////////////////////////////////////////////////////////
  19
  20 template <int DIMENSION> class Meta_Wrapper;
  21
  22 /*********************************************************/
  23 /*                                                       */
  24 /*                    Meta_dTree                         */
  25 /*                                                       */
  26 /*********************************************************/
  27
  28 template <int DIMENSION> class Meta_dTree : public dTree<Wrapper_Noeud<DIMENSION>,Wrapper_Nuage_Noeud<DIMENSION>,DIMENSION>
  29 {
  30 public :
  31         inline int trouve_plus_proche_point(double *node);
  32 };
  33
  34 /*********************************************************/
  35 /*                                                       */
  36 /*                 Meta_Nuage_Maille                     */
  37 /*                                                       */
  38 /*********************************************************/
  39
  40
  41 class Meta_Nuage_Maille : public Wrapper_Nuage_Maille<Wrapper_Med_Connectivity>
  42 {
  43 protected :
  44         Wrapper_Med_Connectivity * connectivite_med;
  45 public :
  46         Meta_Nuage_Maille(CONNECTIVITY * connmed);
  47         Meta_Nuage_Maille():connectivite_med(NULL) {}
  48         ~Meta_Nuage_Maille() {if (connectivite_med) delete connectivite_med;}
  49 };
  50
  51 /*********************************************************/
  52 /*                                                       */
  53 /*                    Meta_Maillage                      */
  54 /*                                                       */
  55 /*********************************************************/
  56
  57
  58 typedef Wrapper_Maillage<Meta_Nuage_Maille> Meta_Maillage;
  59
  60 /*********************************************************/
  61 /*                                                       */
  62 /*                    Meta_Mapping                       */
  63 /*                                                       */
  64 /*********************************************************/
  65
  66 template <int DIMENSION> class Meta_Mapping : public Mapping<Meta_Maillage,Meta_Nuage_Maille,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,DIMENSION>
  67 {
  68 public :
  69         Meta_Mapping(Meta_Wrapper<DIMENSION> * MW):Mapping<Meta_Maillage,Meta_Nuage_Maille,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,DIMENSION>(MW->Get_Maillage(),MW->Get_Nuage_Noeuds(),NULL) {}
  70         Meta_Mapping(Meta_Wrapper<DIMENSION> * MW,Meta_Wrapper<DIMENSION> * TWB):Mapping<Meta_Maillage,Meta_Nuage_Maille,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,DIMENSION>(MW->Get_Maillage(),MW->Get_Nuage_Noeuds(),TWB->Get_Nuage_Noeuds()) {}
  71         void Cree_Mapping(Meta_Wrapper<DIMENSION> * MWB, int flag_convexe) {Mapping<Meta_Maillage,Meta_Nuage_Maille,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,DIMENSION>::Cree_Mapping(MWB->Get_Nuage_Noeuds(),flag_convexe);}
  72         inline int Trouve_Maille_Contenant_Noeud(double * node,int num_maille, int flag_convexe=0);
  73 };
  74
  75 /*********************************************************/
  76 /*                                                       */
  77 /*                    Meta_Wrapper                       */
  78 /*                                                       */
  79 /*********************************************************/
  80
  81
  82 template <int DIMENSION> class Meta_Wrapper
  83 {
  84 protected :
  85         Wrapper_Nuage_Noeud<DIMENSION> * noeuds   ;
  86         Meta_Nuage_Maille              * mailles  ;
  87         Meta_Maillage                  * maillage ;
  88         Wrapper_MED_Field              * champ    ;
  89
  90         void init( ){noeuds=NULL;mailles=NULL;maillage=NULL;champ=NULL;}
  91 public :
  92         Meta_Wrapper():noeuds(NULL),mailles(NULL),maillage(NULL),champ(NULL){}
  93         ~Meta_Wrapper();
  94         inline void Construit_Wrapper_Nuage_Noeud  (     int nn, double * nodes     );
  95         inline void Construit_Wrapper_Nuage_Maille (     CONNECTIVITY * connmed     );
  96         inline void Construit_Wrapper_Maillage     (            void                );
  97         inline void Construit_Wrapper_Champ        ( const FIELD<double> * medfield );
  98         Meta_Wrapper(int nn,double *nodes,CONNECTIVITY *connmed, int flag_maillage=1);
  99         Meta_Wrapper(int nn,double *nodes);
 100         // defaultly, the connectivity (neighbouhood and so like) is built,
 101         // Set flag_mesh to 0 if you don't want these informations to be built
 102         Meta_Wrapper(int nn,double *nodes,CONNECTIVITY *connmed, const FIELD<double> * c,int flag_mesh=1);
 103         // fonctions d'acces sures
 104         inline Wrapper_Nuage_Noeud<DIMENSION> * Get_Nuage_Noeuds  ( void );
 105         inline Meta_Nuage_Maille              * Get_Nuage_Mailles ( void );
 106         inline Meta_Maillage                  * Get_Maillage      ( void );
 107         inline Wrapper_MED_Field              * Get_Champ         ( void );
 108 };
 109
 110 /*********************************************************/
 111 /*                                                       */
 112 /*           Meta_Calcul_Interpolation_Hybride           */
 113 /*                                                       */
 114 /*********************************************************/
 115
 116 template <int DIMENSION> class Meta_Calcul_Interpolation_Hybride : public Calcul_Hybride<Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>
 117 {
 118 public :
 119         Meta_Calcul_Interpolation_Hybride(Meta_Wrapper<DIMENSION> * MW):Calcul_Hybride<Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>(MW->Get_Nuage_Noeuds(),MW->Get_Nuage_Mailles(),MW->Get_Champ()) {}
 120         Valeur<double> operator() (Wrapper_Noeud<DIMENSION> & node, int num_maille){return Calcul_Hybride<Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>::operator()(node,num_maille);}
 121         Valeur<double> operator() (double * node, int num_maille)
 122                 {
 123                 static Wrapper_Noeud<DIMENSION> tmp;
 124                 tmp.positionne(node);
 125                 return Calcul_Hybride<Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>(tmp,num_maille);
 126                 }
 127 };
 128
 129 /*********************************************************/
 130 /*                                                       */
 131 /*         Meta_Calcul_Interpolation_Hybride_P1          */
 132 /*                                                       */
 133 /*********************************************************/
 134
 135 template <int DIMENSION> class Meta_Calcul_Interpolation_Hybride_P1 : public Calcul_Hybride<Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>
 136 {
 137 public :
 138         Meta_Calcul_Interpolation_Hybride_P1(Meta_Wrapper<DIMENSION> * MW)
 139                 {
 140
 141                 Wrapper_Nuage_Noeud<DIMENSION> * nn = MW->Get_Nuage_Noeuds();
 142                 Meta_Nuage_Maille *              nm = MW->Get_Nuage_Mailles();
 143                 Wrapper_MED_Field *              c  = MW->Get_Champ();
 144
 145                 mailles=nm;
 146
 147                 fonctions[MED_TRIA3  ]=new Calcul_Interpolation_Tria3  <Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>(nn,nm,c);
 148                 fonctions[MED_QUAD4  ]=new Calcul_Interpolation_Quad4  <Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>(nn,nm,c);
 149                 fonctions[MED_TETRA4 ]=new Calcul_Interpolation_Tetra4 <Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>(nn,nm,c);
 150                 fonctions[MED_HEXA8  ]=new Calcul_Interpolation_Hexa8  <Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>(nn,nm,c);
 151                 fonctions[MED_PENTA6 ]=new Calcul_Interpolation_Penta6 <Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>(nn,nm,c);
 152                 fonctions[MED_PYRA5  ]=new Calcul_Interpolation_Pyra5  <Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>(nn,nm,c);
 153                 fonctions[MED_TRIA6  ]=fonctions[MED_TRIA3  ];
 154                 fonctions[MED_QUAD8  ]=fonctions[MED_QUAD4  ];
 155                 fonctions[MED_TETRA10]=fonctions[MED_TETRA4 ];
 156                 fonctions[MED_HEXA20 ]=fonctions[MED_HEXA8  ];
 157                 fonctions[MED_PENTA15]=fonctions[MED_PENTA6 ];
 158                 fonctions[MED_PYRA13 ]=fonctions[MED_PYRA5  ];
 159                 }
 160         Valeur<double> operator() (Wrapper_Noeud<DIMENSION> & node, int num_maille){return Calcul_Hybride<Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>::operator()(node,num_maille);}
 161         Valeur<double> operator() (double * node, int num_maille)
 162                 {
 163                 static Wrapper_Noeud<DIMENSION> tmp;
 164                 tmp.positionne(node);
 165                 return Calcul_Hybride<Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>(tmp,num_maille);
 166                 }
 167 };
 168
 169 /*********************************************************/
 170 /*                                                       */
 171 /*              Meta_Calcul_Interpolation_P0             */
 172 /*                                                       */
 173 /*********************************************************/
 174
 175 template <int DIMENSION> class Meta_Calcul_Interpolation_P0 : public Calcul_Interpolation_P0<Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>
 176 {
 177 public :
 178         Meta_Calcul_Interpolation_P0(Meta_Wrapper<DIMENSION> * MW):Calcul_Interpolation_P0<Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>(MW->Get_Nuage_Noeuds(),MW->Get_Nuage_Mailles(),MW->Get_Champ()) {}
 179         Valeur<double> operator() (Wrapper_Noeud<DIMENSION> & node, int num_maille){return Calcul_Interpolation_P0<Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>::operator()(node,num_maille);}
 180         Valeur<double> operator() (double * node, int num_maille)
 181                 {
 182                 static Wrapper_Noeud<DIMENSION> tmp;
 183                 tmp.positionne(node);
 184                 return Calcul_Interpolation_P0<Wrapper_MED_Field,Valeur<double>,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,Meta_Nuage_Maille>(tmp,num_maille);
 185                 }
 186 };
 187
 188 /*********************************************************/
 189 /*                                                       */
 190 /*                   Meta_Interpolateur                  */
 191 /*                                                       */
 192 /*********************************************************/
 193
 194 template <class FONCTEUR, int DIMENSION> class Meta_Interpolateur
 195 {
 196 protected :
 197         FONCTEUR * fct;
 198         Meta_Mapping<DIMENSION> * mapping;
 199         Meta_Wrapper<DIMENSION> * fromWrapper;
 200 public :
 201         Meta_Interpolateur():fct(NULL),mapping(NULL),fromWrapper(NULL) {}
 202         Meta_Interpolateur(Meta_Mapping<DIMENSION> * map, Meta_Wrapper<DIMENSION> * mw):mapping(map),fromWrapper(mw),fct(new FONCTEUR(mw)){}
 203         ~Meta_Interpolateur() {if (fct) delete fct;}
 204         Wrapper_MED_Field Perform_Interpolation(Wrapper_Nuage_Noeud<DIMENSION> * toNodes)
 205                 {
 206                 int i;
 207
 208                 int nbr_composantes = fromWrapper->Get_Champ()->Get_Nbr_Composantes();
 209                 int nbr_valeurs     = toNodes->SIZE();
 210
 211                 double * valeurs=new double[nbr_valeurs*nbr_composantes];
 212
 213                 Wrapper_MED_Field resultat(nbr_valeurs,nbr_composantes,valeurs);
 214
 215                 int nlpp,nmc;
 216
 217                 for (i=0;i<nbr_valeurs;i++)
 218                         {
 219                         cout<<"Interpolation du noeud "<<i<<flush;
 220                         nmc = (*mapping)[i];
 221                         cout<<" | mappé dans la maille "<<nmc<<flush;
 222                         cout<<" | coordonnées = "<<flush<<(*toNodes)[i]<<flush;
 223                         if (nmc>=0)
 224                                 {
 225                                 cout<<" | valeurs qui va etre assignée = "<<flush<<(*fct)((*toNodes)[i],nmc)<<flush;
 226                                 resultat[i]=(*fct)((*toNodes)[i],nmc);
 227                                 }
 228                         else
 229                                 {
 230                                 nlpp = mapping->Get_Noeud_Le_Plus_Proche(i);
 231                                 cout<<" | et dont le point le plus proche a pour numéro : "<<nlpp<<flush;
 232                                 cout<<" | valeurs qui va etre assignée = "<<(*fromWrapper->Get_Champ())[nlpp]<<flush;
 233                                 if (nlpp!=UNDEFINED) resultat[i]=(*fromWrapper->Get_Champ())[nlpp];
 234                                 else
 235                                         {
 236                                         cerr<<"Meta_Interpolateur : Le noeud "<<i+1<<" n'a ni maille contenante, ni point le plus proche"<<flush;
 237                                         exit(-1);
 238                                         }
 239                                 }
 240                         cout<<" | => OK ! "<<endl;
 241                         }
 242
 243                 cout<<"Avant return"<<endl;
 244
 245                 return resultat;
 246
 247                 }
 248 };
 249
 250
 251
 252 //////////////////////////////////////////////////////////////////
 253 ///                                                            ///
 254 ///                            CODE                            ///
 255 ///                                                            ///
 256 //////////////////////////////////////////////////////////////////
 257
 258 /*********************************************************/
 259 /*                                                       */
 260 /*                    Meta_dTree                         */
 261 /*                                                       */
 262 /*********************************************************/
 263
 264 template <int DIMENSION> inline int Meta_dTree<DIMENSION>::trouve_plus_proche_point(double *node)
 265         {
 266         static Wrapper_Noeud<DIMENSION> nodetmp;
 267         nodetmp.positionne(node);
 268         return dTree<Wrapper_Noeud<DIMENSION>,Wrapper_Nuage_Noeud<DIMENSION>,DIMENSION>::trouve_plus_proche_point(Wrapper_Noeud<DIMENSION>(nodetmp));
 269         }
 270
 271 /*********************************************************/
 272 /*                                                       */
 273 /*                 Meta_Nuage_Maille                     */
 274 /*                                                       */
 275 /*********************************************************/
 276
 277 Meta_Nuage_Maille::Meta_Nuage_Maille(CONNECTIVITY * conmed):Wrapper_Nuage_Maille<Wrapper_Med_Connectivity>(connectivite_med=new Wrapper_Med_Connectivity(conmed))
 278         {
 279         }
 280
 281 /*********************************************************/
 282 /*                                                       */
 283 /*                    Meta_Mapping                       */
 284 /*                                                       */
 285 /*********************************************************/
 286
 287 template <int DIMENSION> inline int Meta_Mapping<DIMENSION>::Trouve_Maille_Contenant_Noeud(double * node,int num_maille,int flag_convexe)
 288         {
 289         int interdit=num_maille;
 290         int max_loop=100;
 291         int nme=0;
 292         static Wrapper_Noeud<DIMENSION> nodetmp;
 293         nodetmp.positionne(node);
 294         return Mapping<Meta_Maillage,Meta_Nuage_Maille,Wrapper_Nuage_Noeud<DIMENSION>,Wrapper_Noeud<DIMENSION>,DIMENSION>::Trouve_Maille_Contenant_Point_Mth_Co(nodetmp,num_maille,interdit,max_loop,nme,flag_convexe);
 295         }
 296
 297 /*********************************************************/
 298 /*                                                       */
 299 /*                    Meta_Wrapper                       */
 300 /*                                                       */
 301 /*********************************************************/
 302
 303 template <int DIMENSION>        Meta_Wrapper<DIMENSION>::~Meta_Wrapper()
 304         {
 305         if ( noeuds   ) delete  noeuds   ;
 306         if ( mailles  ) delete  mailles  ;
 307         if ( maillage ) delete  maillage ;
 308         if ( champ    ) delete  champ    ;
 309         }
 310 template <int DIMENSION>        inline void Meta_Wrapper<DIMENSION>::Construit_Wrapper_Nuage_Noeud  ( int nn, double * nodes )
 311         {
 312         if (nodes) noeuds=new Wrapper_Nuage_Noeud<DIMENSION>(nn,nodes);
 313         else
 314                 {
 315                 cerr<<"Meta_Wrapper : Nuage MED_FULL_INTERLACE vide passé en argument au constructeur"<<endl;
 316                 exit(-1);
 317                 }
 318         }
 319 template <int DIMENSION>        inline void Meta_Wrapper<DIMENSION>::Construit_Wrapper_Nuage_Maille ( CONNECTIVITY * connmed )
 320         {
 321         if (connmed) mailles=new Meta_Nuage_Maille(connmed);
 322         else
 323                 {
 324                 cerr<<"Meta_Wrapper : CONNECTIVITY vide passée en argument au constructeur"<<endl;
 325                 exit(-1);
 326                 }
 327         }
 328 template <int DIMENSION>        inline void Meta_Wrapper<DIMENSION>::Construit_Wrapper_Maillage     (         void           )
 329         {
 330         if (mailles==NULL)
 331                 {
 332                 cerr<<"Meta_Wrapper : Le nuage de maille n'a pas été initialisé !"<<endl;
 333                 exit(-1);
 334                 }
 335         if (noeuds==NULL)
 336                 {
 337                 cerr<<"Meta_Wrapper : Le nuage de noeuds n'a pas été initialisé !"<<endl;
 338                 exit(-1);
 339                 }
 340         maillage=new Meta_Maillage(mailles,noeuds->SIZE());
 341         }
 342 template <int DIMENSION>        inline void Meta_Wrapper<DIMENSION>::Construit_Wrapper_Champ        ( const FIELD<double> * medfield )
 343         {
 344         if (medfield) champ=new Wrapper_MED_Field(medfield);
 345         else
 346                 {
 347                 cerr<<"Meta_Wrapper : FIELD MED vide passé en argument au constructeur"<<endl;
 348                 exit(-1);
 349                 }
 350         }
 351 template <int DIMENSION>        Meta_Wrapper<DIMENSION>::Meta_Wrapper(int nn,double *nodes,CONNECTIVITY *connmed, int flag_maillage=1)
 352         {
 353         init();
 354         Construit_Wrapper_Nuage_Noeud(nn,nodes);
 355         Construit_Wrapper_Nuage_Maille(connmed);
 356         if (flag_maillage) Construit_Wrapper_Maillage();
 357         }
 358 template <int DIMENSION>        Meta_Wrapper<DIMENSION>::Meta_Wrapper(int nn,double *nodes,CONNECTIVITY *connmed, const FIELD<double> * c,int flag_maillage=1)
 359         {
 360         init();
 361         Construit_Wrapper_Nuage_Noeud(nn,nodes);
 362         Construit_Wrapper_Nuage_Maille(connmed);
 363         if (flag_maillage) Construit_Wrapper_Maillage();
 364         Construit_Wrapper_Champ(c);
 365         }
 366 template <int DIMENSION>        Meta_Wrapper<DIMENSION>::Meta_Wrapper(int nn,double *nodes)
 367         {
 368         init();
 369         Construit_Wrapper_Nuage_Noeud(nn,nodes);
 370         }
 371 template <int DIMENSION>        inline Wrapper_Nuage_Noeud<DIMENSION> * Meta_Wrapper<DIMENSION>::Get_Nuage_Noeuds  ( void )
 372         {
 373         if (noeuds) return noeuds;
 374         else
 375                 {
 376                 cerr<<"Meta_Wrapper : Nuage noeuds demandé alors qu'il n'est pas construit !"<<endl;
 377                 exit(-1);
 378                 }
 379         }
 380 template <int DIMENSION>        inline Meta_Nuage_Maille              * Meta_Wrapper<DIMENSION>::Get_Nuage_Mailles ( void )
 381         {
 382         if (mailles) return mailles  ;
 383         else
 384                 {
 385                 cerr<<"Meta_Wrapper : Nuage mailles demandé alors qu'il n'est pas construit !"<<endl;
 386                 exit(-1);
 387                 }
 388         }
 389 template <int DIMENSION>        inline Meta_Maillage                  * Meta_Wrapper<DIMENSION>::Get_Maillage      ( void )
 390         {
 391         if (maillage) return maillage ;
 392         else
 393                 {
 394                 cerr<<"Meta_Wrapper : Connectivitée maillage demandée alors qu'elle n'est pas construite !"<<endl;
 395                 exit(-1);
 396                 }
 397         }
 398
 399 template <int DIMENSION>        inline Wrapper_MED_Field              * Meta_Wrapper<DIMENSION>::Get_Champ         ( void )
 400         {
 401         if (champ) return champ;
 402         else
 403                 {
 404                 cerr<<"Meta_Wrapper : Champ demandé alors qu'il n'est pas construit !"<<endl;
 405                 exit(-1);
 406                 }
 407         }
 408
 409 #endif