src/INTERPOLATION/MEDMEM_Interpolation.hxx

   1 # ifndef INTERPOLATION_HXX
   2 # define INTERPOLATION_HXX
   3
   4 //template < class T> class FIELD;
   5 //template < int > class Wrapper_Nuage;
   6 //template < int > class Wrapper_Noeud;
   7 //template <class ,class ,int ,int > class dTree;
   8
   9 #include <vector>
  10 #include "utilities.h"
  11 #include "MEDMEM_Exception.hxx"
  12 #include "MEDMEM_define.hxx"
  13
  14 #include "MEDMEM_InterpolationHighLevelObjects.hxx"
  15 #include "MEDMEM_Mesh.hxx"
  16 #include "MEDMEM_Field.hxx"
  17
  18 namespace MEDMEM {
  19 class MESH;
  20
  21 template <int DIMENSION=3> class INTERPOLATION
  22 {
  23 protected:
  24
  25   FIELD<double> *                 _fromField;
  26   FIELD<double> *                 _toField;
  27   MESH   *                        _fromMesh;
  28   MESH   *                        _toMesh;
  29
  30   Meta_Wrapper<DIMENSION>  *      _fromWrapper;
  31   Meta_Wrapper<DIMENSION>  *      _toWrapper;
  32
  33   Meta_Mapping <DIMENSION> *      _mapping;
  34
  35 // only used when multi timestep are interpolated
  36 // but always coherent
  37   int                             _iType;
  38   int                             _isConvexFromMesh;
  39
  40 public :
  41
  42   void init();
  43
  44   // Initialize INTERPOLATION in order to get :
  45   // 1- the node number in the MESH  <fromMesh> which
  46   //     is the nearest from a given one ( use method : getNearestNode( double * node ) );
  47   // 2- the cell number (if exists) in the MESH <fromMesh> which
  48   //     contains a specified node ( use method : getContainingCell ( double * node) )
  49   INTERPOLATION(const MESH & fromMesh );
  50   // Initialize INTERPOLATION in order to get :
  51   // 1- the complete mapping ( use methode : getMapping() )
  52   // 2- the functionalities above
  53   INTERPOLATION(const MESH & fromMesh,const MESH & toMesh );
  54   // Initialize INTERPOLATION in order to get the interpolation of <field> on <toMesh>
  55   // Moreover, all the others functionalities are so available
  56   INTERPOLATION(const FIELD<double> & fromField, const MESH & toMesh);
  57
  58   ~INTERPOLATION( );
  59
  60   //  Get the node number in the MESH  <fromMesh> which is the nearest from a given one
  61   int getNearestNode    ( double * node );
  62   //  Get the cell number (if exists) in the MESH <fromMesh> which contains a specified node
  63   int getContainingCell ( double * node , int beginingCell=0, int flagIsConvexMesh=0 );
  64   // Get the complete mapping, defaultly, fromMesh is supposed to be non-convex, if it is false, set flagIsConvexMesh to 1
  65   vector<int> getMapping ( int flagIsConvexMesh=0 );
  66   // Get the interpolated field toField
  67   FIELD<double> * interpolate( /*med_interpolation_type*/ int itype,int flagIsConvexFromMesh=0);
  68   // reset the <from> parameters in order not to redo the mapping (if the mesh are identical)
  69   // and then get the interpoated field toField
  70   // this method is specifictly used on multi-timestep (or order number) fields
  71   // it has only to be used after the first step, the interpolation paramaters are the same for every step
  72   FIELD<double> * interpolateNextStep(const FIELD <double> &nextFromField ,int & flagNewMapping);
  73
  74 };
  75
  76 template <int DIMENSION> void INTERPOLATION<DIMENSION>::init() {
  77
  78   const char * LOC = "INTERPOLATION::init(): ";
  79
  80   BEGIN_OF(LOC);
  81   _fromField   = ( FIELD<double> * )           NULL;
  82   _toField     = ( FIELD<double> * )           NULL;
  83   _fromMesh    = ( MESH * )                    NULL;
  84   _toMesh      = ( MESH * )                    NULL;
  85   _fromWrapper = ( Meta_Wrapper<DIMENSION> * ) NULL;
  86   _toWrapper   = ( Meta_Wrapper<DIMENSION> * ) NULL;
  87   _mapping     = ( Meta_Mapping<DIMENSION> * ) NULL;
  88   _iType            = UNDEFINED ;
  89   _isConvexFromMesh = UNDEFINED ;
  90   END_OF(LOC);
  91 }
  92
  93
  94 template <int DIMENSION> INTERPOLATION<DIMENSION>::INTERPOLATION(const MESH & fromMesh ) {
  95
  96   const char * LOC = "INTERPOLATION::INTERPOLATION(MESH * fromMesh ) : ";
  97   BEGIN_OF(LOC);
  98
  99   init();
 100
 101   _fromMesh=const_cast<MESH * > (&fromMesh);
 102
 103   if (! _fromMesh ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"fromMesh is a NULL pointer  !")) ;
 104
 105   int spaceDimension = _fromMesh->getSpaceDimension();
 106   if (spaceDimension != DIMENSION ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The spaceDimension of mesh |" << _fromMesh->getName() << "| is |" << spaceDimension << "| and should be |" << DIMENSION << "|" << endl)) ;
 107
 108   _fromWrapper = new Meta_Wrapper<DIMENSION>(_fromMesh->getNumberOfNodes(),
 109                                              const_cast<double *> (_fromMesh->getCoordinates(MED_FULL_INTERLACE)),
 110                                              const_cast<CONNECTIVITY *> (_fromMesh->getConnectivityptr())
 111                                              );
 112
 113   _mapping     = new  Meta_Mapping<DIMENSION> (_fromWrapper);
 114
 115   END_OF(LOC);
 116 };
 117
 118 template <int DIMENSION> INTERPOLATION<DIMENSION>::INTERPOLATION(const MESH & fromMesh,const MESH & toMesh ) {
 119
 120   const char * LOC = "INTERPOLATION::INTERPOLATION(MESH * fromMesh,,const MESH & toMesh) : ";
 121   BEGIN_OF(LOC);
 122
 123   init();
 124
 125   _fromMesh = const_cast<MESH * > ( &fromMesh );
 126   _toMesh   = const_cast<MESH * > ( &toMesh   );
 127
 128   if (! _fromMesh ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"fromMesh is a NULL pointer  !")) ;
 129   if (! _toMesh   ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toMesh   is a NULL pointer  !")) ;
 130
 131   int fromSpaceDimension = _fromMesh->getSpaceDimension();
 132   int toSpaceDimension   =   _toMesh->getSpaceDimension();
 133
 134   if (fromSpaceDimension != DIMENSION ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The spaceDimension of mesh |" << _fromMesh->getName() << "| is |" << spaceDimension << "| and should be |" << DIMENSION << "|" << endl)) ;
 135   if (  toSpaceDimension != DIMENSION ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The spaceDimension of mesh |" <<   _toMesh->getName() << "| is |" << spaceDimension << "| and should be |" << DIMENSION << "|" << endl)) ;
 136
 137   _fromWrapper = new Meta_Wrapper<DIMENSION>(_fromMesh->getNumberOfNodes(),
 138                                              const_cast<double *> (_fromMesh->getCoordinates(MED_FULL_INTERLACE)),
 139                                              const_cast<CONNECTIVITY *> (_fromMesh->getConnectivityptr())
 140                                              );
 141
 142   _toWrapper   = new Meta_Wrapper<DIMENSION>(_toMesh->getNumberOfNodes(),
 143                                              const_cast<double *> (_toMesh->getCoordinates(MED_FULL_INTERLACE))
 144                                              );
 145
 146   _mapping     = new  Meta_Mapping<DIMENSION> (_fromWrapper);
 147
 148   END_OF(LOC);
 149 };
 150
 151 template <int DIMENSION> INTERPOLATION<DIMENSION>::INTERPOLATION(const FIELD<double> & fromField,const MESH & toMesh) {
 152
 153   const char * LOC = "INTERPOLATION(const FIELD<double> & field,const MESH & toMesh) : ";
 154   BEGIN_OF(LOC);
 155
 156   init();
 157
 158   _toMesh    = const_cast<MESH *>(&toMesh);
 159   _fromField = const_cast<FIELD<double> *>(&fromField);
 160
 161   if ( ! _toMesh    )  throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toMesh is a NULL pointer  !")) ;
 162   if ( ! _fromField )  throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"field is a NULL pointer  !")) ;
 163
 164   _fromMesh = _fromField->getSupport()->getMesh();
 165
 166   if ( ! _fromMesh    )  throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"fromMesh is a NULL pointer  !")) ;
 167
 168   _fromWrapper = new Meta_Wrapper<DIMENSION>(_fromMesh->getNumberOfNodes(),
 169                                              const_cast<double *> (_fromMesh->getCoordinates(MED_EN::MED_FULL_INTERLACE)),
 170                                              const_cast<MEDMEM::CONNECTIVITY *> (_fromMesh->getConnectivityptr()),
 171                                              const_cast<MEDMEM::FIELD<double> *>(_fromField)
 172                                              );
 173
 174
 175   _toWrapper   = new Meta_Wrapper<DIMENSION>(_toMesh->getNumberOfNodes(),
 176                                              const_cast<double *> (_toMesh->getCoordinates(MED_EN::MED_FULL_INTERLACE))
 177                                              );
 178
 179
 180   _mapping     = new  Meta_Mapping<DIMENSION> (_fromWrapper);
 181
 182
 183   END_OF(LOC);
 184 };
 185
 186 template <int DIMENSION> INTERPOLATION<DIMENSION>::~INTERPOLATION()
 187 {
 188   if ( _fromWrapper  ) delete _fromWrapper ;
 189   if ( _toWrapper    ) delete _toWrapper   ;
 190   if ( _mapping      ) delete _mapping     ;
 191 };
 192
 193 template <int DIMENSION> int INTERPOLATION<DIMENSION>::getNearestNode(  double * node ) {
 194
 195   const char * LOC = "INTERPOLATION::getNearestNode( double * node ) ";
 196
 197   BEGIN_OF(LOC);
 198
 199   if ( ! _mapping ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"mapping is a NULL pointer  !")) ;
 200
 201   return _mapping->Donne_dTree()->trouve_plus_proche_point(Wrapper_Noeud<DIMENSION > (node) );
 202
 203   END_OF(LOC);
 204
 205 };
 206
 207 template <int DIMENSION> int INTERPOLATION<DIMENSION>::getContainingCell ( double * node , int beginingCell, int flagIsConvexMesh ) {
 208
 209   const char * LOC = "INTERPOLATION::getContainingCell( double * node ) ";
 210
 211   BEGIN_OF(LOC);
 212
 213   if ( ! _mapping ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"mapping is a NULL pointer  !")) ;
 214
 215   _isConvexFromMesh=flagIsConvexMesh;
 216
 217   _mapping->Cree_Mapping(_toWrapper,_isConvexFromMesh);
 218
 219   return _mapping->Trouve_Maille_Contenant_Noeud(node,beginingCell,flagIsConvexMesh);
 220
 221   END_OF(LOC);
 222
 223 };
 224
 225 template <int DIMENSION> vector<int> INTERPOLATION<DIMENSION>::getMapping ( int flagIsConvexMesh ) {
 226
 227   const char * LOC = "INTERPOLATION::getMapping( ) ";
 228
 229   BEGIN_OF(LOC);
 230
 231   if ( ! _mapping   ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"mapping is a NULL pointer  !")) ;
 232   if ( ! _toWrapper ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toWrapper  is a NULL pointer  !")) ;
 233
 234   _isConvexFromMesh=flagIsConvexMesh;
 235
 236   _mapping->Cree_Mapping(_toWrapper,_isConvexFromMesh);
 237
 238   return _mapping->Get_Mapping();
 239
 240   END_OF(LOC);
 241
 242 };
 243
 244 template <int DIMENSION> FIELD<double> * INTERPOLATION<DIMENSION>::interpolate(int itype,int flagIsConvexFromMesh) {
 245
 246   const char * LOC = "INTERPOLATION::interpolate(int itype,int flagIsConvexFromMesh) ";
 247
 248   BEGIN_OF(LOC);
 249
 250   if ( ! _mapping      ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"mapping is a NULL pointer  !")) ;
 251   if ( ! _toWrapper    ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toWrapper  is a NULL pointer  !")) ;
 252   if ( ! _fromField    ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"fromField  is a NULL pointer  !")) ;
 253
 254   _iType=itype;
 255
 256   _isConvexFromMesh=flagIsConvexFromMesh;
 257
 258   _mapping->Cree_Mapping(_toWrapper,_isConvexFromMesh);
 259
 260   Wrapper_Nuage_Noeud<DIMENSION> * toNodes = _toWrapper->Get_Nuage_Noeuds();
 261
 262   Wrapper_MED_Field resultat;
 263
 264   /*
 265   cout<<"Mapping"<<endl;
 266   _mapping->affiche();
 267   cout<<"Mailles"<<endl;
 268   _fromWrapper->Get_Maillage()->DONNE_POINTEUR_NUAGEMAILLE()->affiche();
 269   cout<<"Noeuds"<<endl;
 270   _fromWrapper->Get_Nuage_Noeuds()->affiche();
 271   */
 272
 273   switch (_iType)
 274     {
 275     case 0 : // INTERPOLATION P0
 276       cout<<"Avant ="<<endl;
 277       resultat=Meta_Interpolateur< Meta_Calcul_Interpolation_P0<DIMENSION>,DIMENSION >(_mapping,_fromWrapper).Perform_Interpolation(toNodes);
 278       break;
 279     case 1 : // INTERPOLATION P-Hybride (Interpole avec la fonction d'interpolation naturelle de la maille contenant le point)
 280       resultat=Meta_Interpolateur< Meta_Calcul_Interpolation_Hybride<DIMENSION>,DIMENSION >(_mapping,_fromWrapper).Perform_Interpolation(toNodes);
 281       break;
 282     case 2 : // INTERPOLATION (P/Q) 1 forcée (Interpole avec la fonction élément fini de la maille de degré 1 -meme si la maille est de degré supérieur-)
 283       resultat=Meta_Interpolateur< Meta_Calcul_Interpolation_Hybride_P1<DIMENSION>, DIMENSION >(_mapping,_fromWrapper).Perform_Interpolation(toNodes);
 284       break;
 285     default :
 286       throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Interpolation type "<<itype<<" not yet implemented !")) ;
 287     }
 288
 289   cout<<"On a fini l'interpolation"<<endl;
 290
 291   _toField = new FIELD<double>;
 292
 293   _toField->setName                   ( _fromField->getName()                   );
 294   _toField->setDescription            ( _fromField->getDescription()            );
 295   _toField->setNumberOfComponents     ( _fromField->getNumberOfComponents()     );
 296   _toField->setNumberOfValues         ( _toMesh   ->getNumberOfNodes()          );
 297   _toField->setComponentsNames        ( _fromField->getComponentsNames()        );
 298   _toField->setComponentsDescriptions ( _fromField->getComponentsDescriptions() );
 299   _toField->setMEDComponentsUnits     ( _fromField->getMEDComponentsUnits()     );
 300   _toField->setIterationNumber        ( _fromField->getIterationNumber()        );
 301   _toField->setTime                   ( _fromField->getTime()                   );
 302   _toField->setOrderNumber            ( _fromField->getOrderNumber()            );
 303   _toField->setValueType              ( MED_EN::MED_REEL64                      );
 304
 305   SUPPORT * mySupport(new SUPPORT(_toMesh,"support",MED_EN::MED_NODE));
 306   _toField->setSupport(mySupport);
 307
 308   _toField->allocValue(_toField->getNumberOfComponents(),_toField->getNumberOfValues());
 309
 310   _toField->setValue( MED_EN::MED_FULL_INTERLACE,resultat.Get_Valeurs());
 311
 312   _toWrapper->Construit_Wrapper_Champ(_toField);
 313
 314   return _toField;
 315
 316   END_OF(LOC);
 317
 318 };
 319
 320 template <int DIMENSION> FIELD<double> * INTERPOLATION<DIMENSION>::interpolateNextStep(const FIELD <double> & nextFromField, int & flagNewMapping) {
 321
 322   const char * LOC = "INTERPOLATION::interpolateNextStep(int itype,int flagIsConvexFromMesh) ";
 323
 324   BEGIN_OF(LOC);
 325
 326   if ( ! _mapping      ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"mapping is a NULL pointer  !"        ));
 327   if ( ! _toWrapper    ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toWrapper  is a NULL pointer  !"     )) ;
 328   if ( ! _fromWrapper  ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"fromWrapper  is a NULL pointer  !"   )) ;
 329   if ( ! _fromField    ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"fromField  is a NULL pointer  !"     )) ;
 330
 331
 332   if ( ! _toField                   ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toField  is a NULL pointer, wrong use of interpolateNextStep"       )) ;
 333   if ( _iType==UNDEFINED            ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"_iType is not defined, wrong use of interpolateNextStep"            )) ;
 334   if ( _isConvexFromMesh==UNDEFINED ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"_isConvexFromMesh is not defined, wrong use of interpolateNextStep" )) ;
 335
 336   // delete _toField; ????????????????????????????
 337
 338   // if the mesh are identical, the mapping is the same, if not, the mapping has to be re-calculated
 339   if (nextFromField.getSupport()->getMesh()->getName()!=_fromMesh->getName())
 340         {
 341
 342         flagNewMapping=1;
 343
 344         delete _mapping;
 345         delete _fromWrapper;
 346
 347         _fromField = const_cast<FIELD<double> *>(&nextFromField);
 348
 349         _fromMesh = _fromField->getSupport()->getMesh();
 350
 351         _fromWrapper = new Meta_Wrapper<DIMENSION>(_fromMesh->getNumberOfNodes(),
 352                                                    const_cast<double *> (_fromMesh->getCoordinates(MED_EN::MED_FULL_INTERLACE)),
 353                                                    const_cast<MEDMEM::CONNECTIVITY *> (_fromMesh->getConnectivityptr()),
 354                                                    const_cast<MEDMEM::FIELD<double> *>(_fromField)
 355                                                    );
 356         _mapping     = new  Meta_Mapping<DIMENSION> (_fromWrapper);
 357
 358         _mapping->Cree_Mapping(_toWrapper,_isConvexFromMesh);
 359
 360         }
 361   else
 362         {
 363
 364         flagNewMapping=0;
 365
 366         _fromField = const_cast<MEDMEM::FIELD<double> *>(&nextFromField);
 367         _fromWrapper->Change_Champ(const_cast<MEDMEM::FIELD<double> *>(_fromField));
 368         }
 369
 370   return interpolate(_iType,_isConvexFromMesh);
 371
 372 };
 373
 374 };
 375
 376 #endif
 377
 378