src/INTERPOLATION/MEDMEM_Interpolation.hxx

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