src/INTERPOLATION/test_MEDMEM_InterpolationFlipBack.cxx

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