src/MEDMEM/INTERPOLATION/MEDMEM_Interpolation.hxx

   1 # ifndef INTERPOLATION_HXX
   2 # define INTERPOLATION_HXX
   3
   4 class MESH;
   5 //template < class T> class FIELD;
   6 //template < int > class Wrapper_Nuage;
   7 //template < int > class Wrapper_Noeud;
   8 //template <class ,class ,int ,int > class dTree;
   9
  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 template <int DIMENSION=3> class INTERPOLATION
  19 {
  20 protected:
  21
  22   FIELD<double> *                 _fromField;
  23   FIELD<double> *                 _toField;
  24   MESH   *                        _fromMesh;
  25   MESH   *                        _toMesh;
  26
  27   Meta_Wrapper<DIMENSION>  *      _fromWrapper;
  28   Meta_Wrapper<DIMENSION>  *      _toWrapper;
  29
  30   Meta_Mapping <DIMENSION> *      _mapping;
  31
  32 public :
  33
  34   void init();
  35
  36   // Initialize INTERPOLATION in order to get :
  37   // 1- the node number in the MESH  <fromMesh> which
  38   //     is the nearest from a given one ( use method : getNearestNode( double * node ) );
  39   // 2- the cell number (if exists) in the MESH <fromMesh> which
  40   //     contains a specified node ( use method : getContainingCell ( double * node) )
  41   INTERPOLATION(const MESH & fromMesh );
  42   // Initialize INTERPOLATION in order to get :
  43   // 1- the complete mapping ( use methode : getMapping() )
  44   // 2- the functionalities above
  45   INTERPOLATION(const MESH & fromMesh,const MESH & toMesh );
  46   // Initialize INTERPOLATION in order to get the interpolation of <field> on <toMesh>
  47   // Moreover, all the others functionalities are so available
  48   INTERPOLATION(const FIELD<double> & fromField, const MESH & toMesh);
  49
  50   ~INTERPOLATION( );
  51
  52   //  Get the node number in the MESH  <fromMesh> which is the nearest from a given one
  53   int getNearestNode    ( double * node );
  54   //  Get the cell number (if exists) in the MESH <fromMesh> which contains a specified node
  55   int getContainingCell ( double * node , int beginingCell=0, int flagIsConvexMesh=0 );
  56   // Get the complete mapping, defaultly, fromMesh is supposed to be non-convex, if it is false, set flagIsConvexMesh to 1
  57   vector<int> getMapping ( int flagIsConvexMesh=0 );
  58   // Get the interpolated field toField
  59   FIELD<double> * interpolate( /*med_interpolation_type*/ int itype,int flagIsConvexFromMesh=0);
  60
  61 };
  62
  63
  64 template <int DIMENSION> void INTERPOLATION<DIMENSION>::init() {
  65
  66   const char * LOC = "INTERPOLATION::init(): ";
  67
  68   BEGIN_OF(LOC);
  69   _fromField   = ( FIELD<double> * )           NULL;
  70   _toField     = ( FIELD<double> * )           NULL;
  71   _fromMesh    = ( MESH * )                    NULL;
  72   _toMesh      = ( MESH * )                    NULL;
  73   _fromWrapper = ( Meta_Wrapper<DIMENSION> * ) NULL;
  74   _toWrapper   = ( Meta_Wrapper<DIMENSION> * ) NULL;
  75   _mapping     = ( Meta_Mapping<DIMENSION> * ) NULL;
  76   END_OF(LOC);
  77 }
  78
  79
  80 template <int DIMENSION> INTERPOLATION<DIMENSION>::INTERPOLATION(const MESH & fromMesh ) {
  81
  82   const char * LOC = "INTERPOLATION::INTERPOLATION(MESH * fromMesh ) : ";
  83   BEGIN_OF(LOC);
  84
  85   init();
  86
  87   _fromMesh=const_cast<MESH * > (&fromMesh);
  88
  89   if (! _fromMesh ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"fromMesh is a NULL pointer  !")) ;
  90
  91   int spaceDimension = _fromMesh->getSpaceDimension();
  92   if (spaceDimension != DIMENSION ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The spaceDimension of mesh |" << _fromMesh->getName() << "| is |" << spaceDimension << "| and should be |" << DIMENSION << "|" << endl)) ;
  93
  94   _fromWrapper = new Meta_Wrapper<DIMENSION>(_fromMesh->getNumberOfNodes(),
  95                                              const_cast<double *> (_fromMesh->getCoordinates(MED_FULL_INTERLACE)),
  96                                              const_cast<CONNECTIVITY *> (_fromMesh->getConnectivityptr())
  97                                              );
  98
  99   _mapping     = new  Meta_Mapping<DIMENSION> (_fromWrapper);
 100
 101   END_OF(LOC);
 102 };
 103
 104 template <int DIMENSION> INTERPOLATION<DIMENSION>::INTERPOLATION(const MESH & fromMesh,const MESH & toMesh ) {
 105
 106   const char * LOC = "INTERPOLATION::INTERPOLATION(MESH * fromMesh,,const MESH & toMesh) : ";
 107   BEGIN_OF(LOC);
 108
 109   init();
 110
 111   _fromMesh = const_cast<MESH * > ( &fromMesh );
 112   _toMesh   = const_cast<MESH * > ( &toMesh   );
 113
 114   if (! _fromMesh ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"fromMesh is a NULL pointer  !")) ;
 115   if (! _toMesh   ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toMesh   is a NULL pointer  !")) ;
 116
 117   int fromSpaceDimension = _fromMesh->getSpaceDimension();
 118   int toSpaceDimension   =   _toMesh->getSpaceDimension();
 119
 120   if (fromSpaceDimension != DIMENSION ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The spaceDimension of mesh |" << _fromMesh->getName() << "| is |" << spaceDimension << "| and should be |" << DIMENSION << "|" << endl)) ;
 121   if (  toSpaceDimension != DIMENSION ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"The spaceDimension of mesh |" <<   _toMesh->getName() << "| is |" << spaceDimension << "| and should be |" << DIMENSION << "|" << endl)) ;
 122
 123   _fromWrapper = new Meta_Wrapper<DIMENSION>(_fromMesh->getNumberOfNodes(),
 124                                              const_cast<double *> (_fromMesh->getCoordinates(MED_FULL_INTERLACE)),
 125                                              const_cast<CONNECTIVITY *> (_fromMesh->getConnectivityptr())
 126                                              );
 127
 128   _toWrapper   = new Meta_Wrapper<DIMENSION>(_toMesh->getNumberOfNodes(),
 129                                              const_cast<double *> (_toMesh->getCoordinates(MED_FULL_INTERLACE))
 130                                              );
 131
 132   _mapping     = new  Meta_Mapping<DIMENSION> (_fromWrapper);
 133
 134   END_OF(LOC);
 135 };
 136
 137 template <int DIMENSION> INTERPOLATION<DIMENSION>::INTERPOLATION(const FIELD<double> & fromField,const MESH & toMesh) {
 138
 139   const char * LOC = "INTERPOLATION::INTERPOLATION(const FIELD<double> & field,const MESH & toMesh) : ";
 140   BEGIN_OF(LOC);
 141
 142   init();
 143
 144   _toMesh    = const_cast<MESH *>(&toMesh);
 145   _fromField = const_cast<FIELD<double> *>(&fromField);
 146
 147   if ( ! _toMesh    )  throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toMesh is a NULL pointer  !")) ;
 148   if ( ! _fromField )  throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"field is a NULL pointer  !")) ;
 149
 150   _fromMesh = _fromField->getSupport()->getMesh();
 151
 152   _fromWrapper = new Meta_Wrapper<DIMENSION>(_fromMesh->getNumberOfNodes(),
 153                                              const_cast<double *> (_fromMesh->getCoordinates(MED_FULL_INTERLACE)),
 154                                              const_cast<CONNECTIVITY *> (_fromMesh->getConnectivityptr()),
 155                                              const_cast<FIELD<double> *>(_fromField)
 156                                              );
 157
 158   _toWrapper   = new Meta_Wrapper<DIMENSION>(_toMesh->getNumberOfNodes(),
 159                                              const_cast<double *> (_toMesh->getCoordinates(MED_FULL_INTERLACE))
 160                                              );
 161
 162   _mapping     = new  Meta_Mapping<DIMENSION> (_fromWrapper);
 163
 164   END_OF(LOC);
 165 };
 166
 167 template <int DIMENSION> INTERPOLATION<DIMENSION>::~INTERPOLATION()
 168 {
 169   if ( _fromWrapper  ) delete _fromWrapper ;
 170   if ( _toWrapper    ) delete _toWrapper   ;
 171   if ( _mapping      ) delete _mapping     ;
 172 };
 173
 174 template <int DIMENSION> int INTERPOLATION<DIMENSION>::getNearestNode(  double * node ) {
 175
 176   const char * LOC = "INTERPOLATION::getNearestNode( double * node ) ";
 177
 178   BEGIN_OF(LOC);
 179
 180   if ( ! _mapping ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"mapping is a NULL pointer  !")) ;
 181
 182   return _mapping->Donne_dTree()->trouve_plus_proche_point(Wrapper_Noeud<DIMENSION > (node) );
 183
 184   END_OF(LOC);
 185
 186 };
 187
 188 template <int DIMENSION> int INTERPOLATION<DIMENSION>::getContainingCell ( double * node , int beginingCell, int flagIsConvexMesh ) {
 189
 190   const char * LOC = "INTERPOLATION::getContainingCell( double * node ) ";
 191
 192   BEGIN_OF(LOC);
 193
 194   if ( ! _mapping ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"mapping is a NULL pointer  !")) ;
 195
 196   return _mapping->Trouve_Maille_Contenant_Noeud(node,beginingCell,flagIsConvexMesh);
 197
 198   END_OF(LOC);
 199
 200 };
 201
 202 template <int DIMENSION> vector<int> INTERPOLATION<DIMENSION>::getMapping ( int flagIsConvexMesh ) {
 203
 204   const char * LOC = "INTERPOLATION::getMapping( ) ";
 205
 206   BEGIN_OF(LOC);
 207
 208   if ( ! _mapping   ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"mapping is a NULL pointer  !")) ;
 209   if ( ! _toWrapper ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toWrapper  is a NULL pointer  !")) ;
 210
 211   mapping_->Cree_Mapping(_toWrapper,flag_convexe);
 212
 213   return _mapping->Get_Mapping();
 214
 215   END_OF(LOC);
 216
 217 };
 218
 219 template <int DIMENSION> FIELD<double> * INTERPOLATION<DIMENSION>::interpolate( /*med_interpolation_type*/ int itype,int flagIsConvexFromMesh) {
 220
 221   const char * LOC = "INTERPOLATION::interpolate( /*med_interpolation_type*/ int itype) ";
 222
 223   BEGIN_OF(LOC);
 224
 225   if ( ! _mapping   ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"mapping is a NULL pointer  !")) ;
 226   if ( ! _toWrapper ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"toWrapper  is a NULL pointer  !")) ;
 227   if ( ! _fromField ) throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"fromField  is a NULL pointer  !")) ;
 228
 229   _mapping->Cree_Mapping(_toWrapper,flagIsConvexFromMesh);
 230
 231   Wrapper_Nuage_Noeud<DIMENSION> * toNodes = _toWrapper->Get_Nuage_Noeuds();
 232
 233   Wrapper_MED_Field resultat;
 234
 235 cout<<"On commence l'interpolation"<<endl;
 236
 237   switch (itype)
 238         {
 239         case 0 : // INTERPOLATION P0
 240                 cout<<"Avant ="<<endl;
 241                 resultat=Meta_Interpolateur< Meta_Calcul_Interpolation_P0<DIMENSION>,DIMENSION >(_mapping,_fromWrapper).Perform_Interpolation(toNodes);
 242                 break;
 243         case 1 : // INTERPOLATION P-Hybride (Interpole avec la fonction d'interpolation naturelle de la maille contenant le point)
 244                 resultat=Meta_Interpolateur< Meta_Calcul_Interpolation_Hybride<DIMENSION>,DIMENSION >(_mapping,_fromWrapper).Perform_Interpolation(toNodes);
 245                 break;
 246         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-)
 247                 resultat=Meta_Interpolateur< Meta_Calcul_Interpolation_Hybride_P1<DIMENSION>, DIMENSION >(_mapping,_fromWrapper).Perform_Interpolation(toNodes);
 248                 break;
 249         default :
 250                 throw MEDEXCEPTION(LOCALIZED(STRING(LOC)<<"Interpolation type "<<itype<<" not yet implemented !")) ;
 251         }
 252
 253 cout<<"On a fini l'interpolation"<<endl;
 254
 255   _toField = new FIELD<double>;
 256
 257   _toField->setName                   ( _fromField->getName()+"XXX"             );
 258   _toField->setDescription            ( _fromField->getDescription()            );
 259   _toField->setNumberOfComponents     ( _fromField->getNumberOfComponents()     );
 260   _toField->setNumberOfValues         ( _toMesh   ->getNumberOfNodes()          );
 261   _toField->setComponentsNames        ( _fromField->getComponentsNames()        );
 262   _toField->setComponentsDescriptions ( _fromField->getComponentsDescriptions() );
 263   _toField->setMEDComponentsUnits     ( _fromField->getMEDComponentsUnits()     );
 264   _toField->setIterationNumber        ( _fromField->getIterationNumber()        );
 265   _toField->setTime                   ( _fromField->getTime()                   );
 266   _toField->setOrderNumber            ( _fromField->getOrderNumber()            );
 267   _toField->setValueType              ( MED_EN::MED_REEL64                      );
 268
 269   SUPPORT * mySupport(new SUPPORT(_toMesh,"support",MED_NODE));
 270   _toField->setSupport(mySupport);
 271
 272   _toField->allocValue(_toField->getNumberOfComponents(),_toField->getNumberOfValues());
 273
 274   _toField->setValue(MED_FULL_INTERLACE,resultat.Get_Valeurs());
 275
 276   _toWrapper->Construit_Wrapper_Champ(_toField);
 277
 278   return _toField;
 279
 280   END_OF(LOC);
 281
 282 };
 283
 284 #endif
 285
 286