src/INTERPOLATION/test_MEDMEM_InterpolationFlipBack.cxx

   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 #include "MEDMEM_Exception.hxx"
  23 #include "MEDMEM_define.hxx"
  24
  25 #include "MEDMEM_Med.hxx"
  26 #include "MEDMEM_Field.hxx"
  27 #include "MEDMEM_Mesh.hxx"
  28 #include "MEDMEM_Interpolation.hxx"
  29
  30 #include <deque>
  31
  32 #include "stdio.h"
  33
  34 using namespace MEDMEM;
  35
  36 // pour gestion timings
  37 #include "time.h"
  38
  39 #define RUN(procedure) {double t0,t1;cout<<"# =============> TEMPS D'EXECUTION A PARTIR D'ICI "<<endl<<#procedure<<endl;t0=CPUtime();procedure;t1=CPUtime();cout<<"# ================> TEMPS D'EXECUTION : "<<t1-t0<<endl;}
  40 #define TIMORIZE(procedure,t) {double t0,t1;t0=CPUtime();procedure;t1=CPUtime();t=t1-t0;}
  41
  42 double CPUtime()
  43         {
  44         #ifdef SYSTIMES
  45         struct tms buf;
  46         if (times(&buf)!=-1)
  47                 return ((double)buf.tms_utime+(double)buf.tms_stime)/(long) sysconf(_SC_CLK_TCK);
  48         else
  49         #endif
  50         return ((double) clock())/CLOCKS_PER_SEC;
  51         }
  52
  53 double ABS(Valeur<double> v)
  54         {
  55         double tmp=0;
  56         int i;
  57         for (i=0;i<v.SIZE();i++) tmp+=fabs(v[i]);
  58         return tmp;
  59         }
  60
  61 double ABS(Valeur<double> v1,Valeur<double> v2)
  62         {
  63         double tmp=0;
  64         int i;
  65         for (i=0;i<v1.SIZE();i++) tmp+=fabs(v1[i]-v2[i]);
  66         return tmp;
  67         }
  68
  69
  70
  71 // effectue la soustraction du premier moins le second et stocke le résultat dans le premier
  72 void Flipback(FIELD<double> * firstField, FIELD<double> * secondField)
  73         {
  74         Wrapper_MED_Field first  ( firstField);
  75         Wrapper_MED_Field second (secondField);
  76         int nbr_valeurs_first  =  first.Get_Nbr_Valeurs();
  77         int nbr_valeurs_second = second.Get_Nbr_Valeurs();
  78
  79         double max1    = 0;
  80         double max2    = 0;
  81
  82         double min1    = ABS(first[0]);
  83         double min2    = ABS(second[0]);
  84
  85         int imax1,imax2;
  86
  87         double tmp;
  88
  89         int i;
  90
  91         //cout<<first<<endl;
  92         //int tyty;cin>>tyty;
  93
  94         if (nbr_valeurs_first!=nbr_valeurs_second)
  95                 {
  96                 cerr<<"Les champs à soustraire n'ont pas le meme nombre de valeurs"<<endl;
  97                 exit(-1);
  98                 }
  99
 100         imax1=0;
 101         for (i=0;i<nbr_valeurs_first;i++)
 102                 {
 103                 tmp=ABS(first[i]);
 104                 //cout<<"tmp 1 ["<<i<<"] = "<<tmp<<endl;
 105                 if (tmp>max1)
 106                         {
 107                         imax1=i;
 108                         max1=tmp;
 109                         }
 110                 if (tmp<min1) min1=tmp;
 111                 }
 112
 113         imax2=0;
 114         for (i=0;i<nbr_valeurs_first;i++)
 115                 {
 116                 tmp=ABS(second[i]);
 117                 if (tmp>max2)
 118                         {
 119                         imax2=i;
 120                         max2=tmp;
 121                         }
 122                 if (tmp<min2) min2=tmp;
 123                 }
 124
 125         for (i=0;i<nbr_valeurs_first;i++)
 126                 {
 127                 first[i]=ABS(first[i],second[i]);
 128                 }
 129
 130         double maxdiff=ABS(first[0]);
 131         double mindiff=ABS(first[0]);
 132
 133         for (i=0;i<nbr_valeurs_first;i++)
 134                 {
 135                 tmp=ABS(first[i]);
 136                 if (tmp>maxdiff) maxdiff=tmp;
 137                 if (tmp<mindiff) mindiff=tmp;
 138                 }
 139
 140         cout<<endl;
 141         cout<<"/////////////////////////////////////////////////////////////////////////"<<endl;
 142         cout<<"/////// max1    = "<<max1<<endl;
 143         cout<<"/////// min1    = "<<min1<<endl;
 144         cout<<"/////// Maximum First  atteint pour i = "<<imax1<<endl;
 145         cout<<"/////// max2    = "<<max2<<endl;
 146         cout<<"/////// min2    = "<<min2<<endl;
 147         cout<<"/////// Maximum Second atteint pour i = "<<imax2<<endl;
 148         cout<<"/////// maxdiff = "<<maxdiff<<endl;
 149         cout<<"/////// mindiff = "<<mindiff<<endl;
 150         cout<<"/////////////////////////////////////////////////////////////////////////"<<endl;
 151
 152         }
 153
 154
 155 #include "MEDMEM_WrapperCells.hxx"
 156
 157 int main () {
 158   const char * fromFileName   = "ResultatSyrthes.med";
 159   const char * toFileName     = "MaillageAster.med";
 160   const char * resultFileName = "ResultatFlipback.med";
 161
 162   const char * fromFieldName  = "THERDEP_TEMP____________________";
 163
 164   const char * fromMeshName   = "MA";
 165   const char * toMeshName     = "MAILLAGE_IDEAS";
 166   int handle;
 167
 168   try {
 169
 170     string flag="================[MAIN MESSAGE_MEDS]================> ";
 171
 172     cout<<flag<<"Lecture de la structure MED               : "<<flush;
 173     MED             fromMED     (MED_DRIVER,fromFileName);
 174     cout<<"OK !"<<endl;
 175
 176     // Utilisation completement débile, on ne devrait pas avoir a faire l'appel suivant
 177     fromMED.updateSupport();
 178
 179     cout<<flag<<"Lecture du Mailllage Cible                : "<<flush;
 180     MESH            toMesh      (MED_DRIVER,toFileName,toMeshName);
 181     cout<<"OK !"<<endl;
 182
 183     cout<<flag<<"Lecture des pas de temps                  : "<<flush;
 184     deque<DT_IT_> pasDeTemps=fromMED.getFieldIteration (fromFieldName);
 185     cout<<"OK !"<<endl;
 186
 187     deque<DT_IT_>::const_iterator currentStep;
 188
 189     INTERPOLATION<3>  * interFromTo ;
 190     INTERPOLATION<3>  * interToFrom ;
 191     FIELD<double>   * toField   ;
 192     FIELD<double>   * toToField ;
 193     int flagNewMappingFromTo = 0;
 194     int flagNewMappingToFrom = 0;
 195
 196     for (currentStep=pasDeTemps.begin();currentStep!=pasDeTemps.end();currentStep++)
 197         {
 198         cout<<flag<<"Traitement du Step ( "<<flush<<(*currentStep).dt<<" ; "<<(*currentStep).it<<" )  : "<<endl;
 199
 200         cout<<flag<<"Lecture du FIELD_ "<<flush;
 201         FIELD_ * fromField_ = fromMED.getField(fromFieldName,(*currentStep).dt,(*currentStep).it);
 202         cout<<"OK !"<<endl;
 203
 204         cout<<flag<<"Transtypage en FIELD                      : "<<flush;
 205         FIELD<double> * fromField = dynamic_cast<FIELD<double> *>(fromField_);
 206         cout<<"OK !"<<endl;
 207
 208         if (currentStep==pasDeTemps.begin())
 209                 {
 210                 //Utilisation completement débile, on ne devrait pas avoir a faire l'appel suivant
 211                 RUN(fromField->getSupport()->getMesh()->read());
 212                 }
 213
 214         MESH * fromMesh = fromField->getSupport()->getMesh();
 215
 216         cout<<flag<<"Lecture des valeurs du FIELD              : "<<flush;
 217         RUN(fromField->read());
 218         cout<<"OK !"<<endl;
 219
 220         if (currentStep==pasDeTemps.begin())
 221                 {
 222                 cout<<flag<<"Préparation de l'interpolation DIRECTE pour le premier pas de temps  : "<<flush;
 223                 RUN(interFromTo = new INTERPOLATION<3>(*fromField,toMesh));
 224                 cout<<"OK !"<<endl;
 225                 cout<<flag<<"Interpolation effective DIRECTE du premier pas de temps              : "<<flush;
 226                 RUN(toField = interFromTo->interpolate(1,1));
 227                 cout<<"OK !"<<endl;
 228                 cout<<flag<<"Préparation de l'interpolation INVERSE pour le premier pas de temps  : "<<flush;
 229                 RUN(interToFrom = new INTERPOLATION<3>(*toField,*fromMesh));
 230                 cout<<"OK !"<<endl;
 231                 cout<<flag<<"Interpolation effective INVERSE du premier pas de temps              : "<<flush;
 232                 RUN(toToField = interToFrom->interpolate(1,1));
 233                 cout<<"OK !"<<endl;
 234                 }
 235         else
 236                 {
 237                 cout<<flag<<"Interpolation nextStep DIRECTE             : "<<flush;
 238                 RUN(toField = interFromTo->interpolateNextStep(*fromField,flagNewMappingFromTo));
 239                 cout<<"OK !"<<endl;
 240                 cout<<flag<<"Interpolation nextStep INVERSE             : "<<flush;
 241                 RUN(toToField = interToFrom->interpolateNextStep(*toField,flagNewMappingToFrom));
 242                 cout<<"OK !"<<endl;
 243                 }
 244
 245         cout<<flag<<"Calcul du flip back : "<<flush;
 246         Flipback(toToField,fromField);
 247         cout<<"OK !"<<endl;
 248
 249         cout<<flag<<"Creation du driver d'écriture Field       : "<<flush;
 250         toToField->addDriver(MED_DRIVER,resultFileName,toToField->getName());
 251         cout<<"OK !"<<endl;
 252
 253         cout<<flag<<"Ecriture du Field résultat                : "<<flush;
 254         toToField->write();
 255         cout<<"OK !"<<endl;
 256
 257         if (flagNewMappingToFrom==1)
 258                 {
 259                 cout<<flag<<"Creation du driver d'écriture Mesh        : "<<flush;
 260                 handle = fromMesh->addDriver(MED_DRIVER,resultFileName,fromMesh->getName()) ;
 261                 cout<<"OK !"<<endl;
 262
 263                 cout<<flag<<"Ecriture du Mesh résultat                 : "<<flush;
 264                 fromMesh->write(handle);
 265                 cout<<"OK !"<<endl;
 266                 }
 267         }
 268
 269   } catch (MEDEXCEPTION& ex){
 270     MESSAGE_MED(ex.what()) ;
 271   }
 272 }