src/MEDMEM/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
 138   const char * fromFileName   = "ResultatSyrthes.med";
 139   const char * toFileName     = "MaillageAster.med";
 140   const char * resultFileName = "ResultatFlipback.med";
 141
 142   const char * fromFieldName  = "THERDEP_TEMP____________________";
 143
 144   const char * fromMeshName   = "MA";
 145   const char * toMeshName     = "MAILLAGE_IDEAS";
 146   int handle;
 147
 148   try {
 149
 150     string flag="================[MAIN MESSAGES]================> ";
 151
 152     cout<<flag<<"Lecture de la structure MED               : "<<flush;
 153     MED             fromMED     (MED_DRIVER,fromFileName);
 154     cout<<"OK !"<<endl;
 155
 156     // Utilisation completement débile, on ne devrait pas avoir a faire l'appel suivant
 157     fromMED.updateSupport();
 158
 159     cout<<flag<<"Lecture du Mailllage Cible                : "<<flush;
 160     MESH            toMesh      (MED_DRIVER,toFileName,toMeshName);
 161     cout<<"OK !"<<endl;
 162
 163     cout<<flag<<"Lecture des pas de temps                  : "<<flush;
 164     deque<DT_IT_> pasDeTemps=fromMED.getFieldIteration (fromFieldName);
 165     cout<<"OK !"<<endl;
 166
 167     deque<DT_IT_>::const_iterator currentStep;
 168
 169     INTERPOLATION<3>  * interFromTo ;
 170     INTERPOLATION<3>  * interToFrom ;
 171     FIELD<double>   * toField   ;
 172     FIELD<double>   * toToField ;
 173     int flagNewMappingFromTo = 0;
 174     int flagNewMappingToFrom = 0;
 175
 176     for (currentStep=pasDeTemps.begin();currentStep!=pasDeTemps.end();currentStep++)
 177         {
 178         cout<<flag<<"Traitement du Step ( "<<flush<<(*currentStep).dt<<" ; "<<(*currentStep).it<<" )  : "<<endl;
 179
 180         cout<<flag<<"Lecture du FIELD_ "<<flush;
 181         FIELD_ * fromField_ = fromMED.getField(fromFieldName,(*currentStep).dt,(*currentStep).it);
 182         cout<<"OK !"<<endl;
 183
 184         cout<<flag<<"Transtypage en FIELD                      : "<<flush;
 185         FIELD<double> * fromField = dynamic_cast<FIELD<double> *>(fromField_);
 186         cout<<"OK !"<<endl;
 187
 188         if (currentStep==pasDeTemps.begin())
 189                 {
 190                 //Utilisation completement débile, on ne devrait pas avoir a faire l'appel suivant
 191                 RUN(fromField->getSupport()->getMesh()->read());
 192                 }
 193
 194         MESH * fromMesh = fromField->getSupport()->getMesh();
 195
 196         cout<<flag<<"Lecture des valeurs du FIELD              : "<<flush;
 197         RUN(fromField->read());
 198         cout<<"OK !"<<endl;
 199
 200         if (currentStep==pasDeTemps.begin())
 201                 {
 202                 cout<<flag<<"Préparation de l'interpolation DIRECTE pour le premier pas de temps  : "<<flush;
 203                 RUN(interFromTo = new INTERPOLATION<3>(*fromField,toMesh));
 204                 cout<<"OK !"<<endl;
 205                 cout<<flag<<"Interpolation effective DIRECTE du premier pas de temps              : "<<flush;
 206                 RUN(toField = interFromTo->interpolate(1,1));
 207                 cout<<"OK !"<<endl;
 208                 cout<<flag<<"Préparation de l'interpolation INVERSE pour le premier pas de temps  : "<<flush;
 209                 RUN(interToFrom = new INTERPOLATION<3>(*toField,*fromMesh));
 210                 cout<<"OK !"<<endl;
 211                 cout<<flag<<"Interpolation effective INVERSE du premier pas de temps              : "<<flush;
 212                 RUN(toToField = interToFrom->interpolate(1,1));
 213                 cout<<"OK !"<<endl;
 214                 }
 215         else
 216                 {
 217                 cout<<flag<<"Interpolation nextStep DIRECTE             : "<<flush;
 218                 RUN(toField = interFromTo->interpolateNextStep(*fromField,flagNewMappingFromTo));
 219                 cout<<"OK !"<<endl;
 220                 cout<<flag<<"Interpolation nextStep INVERSE             : "<<flush;
 221                 RUN(toToField = interToFrom->interpolateNextStep(*toField,flagNewMappingToFrom));
 222                 cout<<"OK !"<<endl;
 223                 }
 224
 225         cout<<flag<<"Calcul du flip back : "<<flush;
 226         Flipback(toToField,fromField);
 227         cout<<"OK !"<<endl;
 228
 229         cout<<flag<<"Creation du driver d'écriture Field       : "<<flush;
 230         toToField->addDriver(MED_DRIVER,resultFileName,toToField->getName());
 231         cout<<"OK !"<<endl;
 232
 233         cout<<flag<<"Ecriture du Field résultat                : "<<flush;
 234         toToField->write();
 235         cout<<"OK !"<<endl;
 236
 237         if (flagNewMappingToFrom==1)
 238                 {
 239                 cout<<flag<<"Creation du driver d'écriture Mesh        : "<<flush;
 240                 handle = fromMesh->addDriver(MED_DRIVER,resultFileName,fromMesh->getName()) ;
 241                 cout<<"OK !"<<endl;
 242
 243                 cout<<flag<<"Ecriture du Mesh résultat                 : "<<flush;
 244                 fromMesh->write(handle);
 245                 cout<<"OK !"<<endl;
 246                 }
 247         }
 248
 249   } catch (MEDEXCEPTION& ex){
 250     MESSAGE(ex.what()) ;
 251   }
 252 }