src/MEDMEM/test_operation_fielddouble.cxx

   1 // Programme de test des operations sur les champs
   2
   3 #include <string>
   4 #include <iostream>
   5 #include <iomanip>
   6 #include <cmath>
   7
   8 #include "MEDMEM_Exception.hxx"
   9 #include "MEDMEM_Mesh.hxx"
  10 #include "MEDMEM_Family.hxx"
  11 #include "MEDMEM_Group.hxx"
  12
  13 #include "MEDMEM_MedMeshDriver.hxx"
  14 #include "MEDMEM_MedFieldDriver.hxx"
  15 #include "MEDMEM_Support.hxx"
  16 #include "MEDMEM_Field.hxx"
  17 #include "MEDMEM_define.hxx"
  18
  19 double myfunction1(double x)
  20 {
  21     return 0.25*(x-1.0);
  22 }
  23
  24
  25 using namespace std;
  26 void affiche_field_(FIELD_ * myField, const SUPPORT * mySupport)
  27 {
  28   cout << "Field "<< myField->getName() << " : " <<myField->getDescription() <<  endl ;
  29   int NumberOfComponents = myField->getNumberOfComponents() ;
  30   cout << "- Nombre de composantes : "<< NumberOfComponents << endl ;
  31   cout << "- Nombre de valeurs     : "<< myField->getNumberOfValues() << endl ;
  32   for (int i=1; i<NumberOfComponents+1; i++) {
  33     cout << "  - composante "<<i<<" :"<<endl ;
  34     cout << "      - nom         : "<<myField->getComponentName(i)<< endl;
  35     cout << "      - description : "<<myField->getComponentDescription(i) << endl;
  36     cout << "      - units       : "<<myField->getMEDComponentUnit(i) << endl;
  37   }
  38   cout << "- iteration :" << endl ;
  39   cout << "    - numero : " << myField->getIterationNumber()<< endl  ;
  40   cout << "    - ordre  : " << myField->getOrderNumber()<< endl  ;
  41   cout << "    - temps  : " << myField->getTime()<< endl  ;
  42
  43   cout << "- Type : " << myField->getValueType()<< endl;
  44
  45   cout << "- Adresse support : " << mySupport << endl;
  46 }
  47
  48 void affiche_fieldT(FIELD<double> * myField, const SUPPORT * mySupport)
  49 {
  50   affiche_field_((FIELD_ *) myField, mySupport);
  51
  52   cout << "- Valeurs :"<<endl;
  53   int NumberOf = mySupport->getNumberOfElements(MED_ALL_ELEMENTS);
  54   int NumberOfComponents = myField->getNumberOfComponents() ;
  55
  56   for (int i=1; i<NumberOf+1; i++) {
  57     const double * value = myField->getValueI(MED_FULL_INTERLACE,i) ;
  58     for (int j=0; j<NumberOfComponents; j++)
  59       cout << value[j]<< " ";
  60     cout<<endl;
  61   }
  62   cout << endl;
  63   cout << "Norme euclidienne : " << myField->norm2() << endl;
  64   cout << "Norme max         : " << myField->normMax() << endl;
  65   try
  66   {
  67       for (int i=1; i<=myField->getNumberOfComponents(); ++i)
  68             cout << "Norme L2 - comp=" << i << " : " << myField->normL2(i) << endl;
  69       cout << "Norme L2          : " << myField->normL2() << endl;
  70
  71       for (int i=1; i<=myField->getNumberOfComponents(); ++i)
  72             cout << "Norme L1 - comp=" << i << " : " << myField->normL1(i) << endl;
  73       cout << "Norme L1          : " << myField->normL1() << endl;
  74   }
  75   catch (MEDEXCEPTION &ex)
  76   {
  77       cout << ex.what() << endl;
  78   }
  79 }
  80
  81 void affiche_valeur_field(const FIELD<double>& f)
  82 {
  83     const int tailleMax=12;
  84     const int taille=f.getNumberOfValues()*f.getNumberOfComponents();
  85     const double * value=f.getValue(f.getvalue()->getMode());
  86     if(taille<=tailleMax)
  87         for(int i=0;i<taille;i++)
  88             cout << setw(3) << value[i] << " ";
  89     else
  90     {
  91         for(int i=0; i<tailleMax/2; ++i)
  92             cout << setw(3) << value[i] << " ";
  93         cout << "    ...    ";
  94         for(int i=taille-tailleMax/2 ; i<taille; ++i)
  95             cout << setw(3) << value[i] << " ";
  96     }
  97 }
  98
  99 void checkOperation(const FIELD<double>& resOp, const FIELD<double>& f1, const FIELD<double>& f2,
 100         char Op, const char* intitule, int verbose)
 101 {
 102     int res=0;
 103
 104     // get pointers to inside arrays of values
 105     medModeSwitch mode=resOp.getvalue()->getMode();
 106     const double * value=resOp.getValue(mode);
 107     const double * value1=f1.getValue(mode);
 108     const double * value2=f2.getValue(mode);
 109     const int size=f1.getNumberOfValues()*f1.getNumberOfComponents(); // size of field1
 110
 111     // check size compatibility
 112     if(f1.getNumberOfValues()*f1.getNumberOfComponents()!=size ||
 113             resOp.getNumberOfValues()*resOp.getNumberOfComponents()!=size)
 114         res=1;
 115
 116     if(!res)
 117     {
 118         switch(Op)
 119         {
 120             case '+':
 121                 for(int i=0; i!=size; ++i)
 122                     if(value[i]!=value1[i]+value2[i])
 123                         res+=1;
 124                 break;
 125             case '-':
 126                 for(int i=0; i!=size; ++i)
 127                     if(value[i]!=value1[i]-value2[i])
 128                         res+=1;
 129                 break;
 130             case 'n':
 131                 for(int i=0; i!=size; ++i)
 132                     if(value[i]!=-value1[i])
 133                         res+=1;
 134                 break;
 135             case '*':
 136                 for(int i=0; i!=size; ++i)
 137                     if(value[i]!=value1[i]*value2[i])
 138                         res+=1;
 139                 break;
 140             case '/':
 141                 for(int i=0; i!=size; ++i)
 142                     if(value2[i]!=0.0)
 143                         if(value[i]!=value1[i]/value2[i])
 144                             res+=1;
 145                 break;
 146             case '=':
 147                 for(int i=0; i!=size; ++i)
 148                     if(value[i]!=value2[i])
 149                         res+=1;
 150                 break;
 151             case 'a':
 152                 for(int i=0; i!=size; ++i)
 153                     if(value[i]!=value1[i]+value2[i]*value2[i])
 154                         res+=1;
 155                 break;
 156         }
 157
 158     }
 159
 160     if (verbose)
 161         cout << endl << intitule << "[";
 162     cout << res;
 163     if (verbose)
 164     {
 165         cout << "] : ";
 166         affiche_valeur_field(resOp);
 167     }
 168     else
 169         cout << endl;
 170 }
 171
 172 int main (int argc, char ** argv)
 173 {
 174     /* process the arguments */
 175     int verbose=0;  //  verbose=1 if the verbose mode is selected
 176     int res=0; // unit test result
 177     int ntest=0;  // numéro du test
 178
 179     if (argc>=2 && !strcmp(argv[1],"-v"))
 180         verbose=1;
 181
 182     if (argc != 4+verbose)
 183     {
 184         cerr << "Usage : " << argv[0]
 185         << "[-v] filename meshname fieldname" << endl << endl
 186         << "-> tests field's operations on the FIELD<double> fieldname" << endl
 187         << "Use optional option -v to select verbose mode" << endl;
 188         exit(-1);
 189     }
 190     string filename  = argv[verbose+1];
 191     string meshname  = argv[verbose+2];
 192     string fieldname = argv[verbose+3];
 193
 194     /* read MESH, SUPPORT and FIELDS */
 195     MESH * myMesh = new MESH(MED_DRIVER,filename,meshname);
 196     SUPPORT * mySupport;
 197     FIELD<double> * myField1;
 198     try
 199     {
 200         mySupport = new SUPPORT(myMesh,"Support on all Cells",MED_CELL);
 201         myField1 = new FIELD<double>(mySupport,MED_DRIVER,filename,fieldname) ;
 202         myField1->setValueType(MED_REEL64);
 203     }
 204     catch (MEDEXCEPTION &ex)
 205     {
 206         // field wasn't found on cells, try on nodes
 207         delete mySupport ;
 208         mySupport = new SUPPORT(myMesh,"On_all_node",MED_NODE);
 209         try
 210         {
 211             myField1 = new FIELD<double>(mySupport,MED_DRIVER,filename,fieldname) ;
 212         }
 213         catch (...)
 214         {
 215             cout << "Field double " << fieldname << " not found !!!" << endl ;
 216             exit (-1) ;
 217         }
 218     }
 219     FIELD<double> * myField2 = new FIELD<double>(* myField1);
 220     FIELD<double> myFieldPlus = *myField1 + *myField2;
 221     if(verbose)
 222     {
 223         // affichage des nprmes,des champs f1, f2, scalarProduct(f1,f2) et f1+f2
 224         FIELD<double>* myField1_vol=myField1->getSupport()->getMesh()->getVolume(myField1->getSupport());
 225         cout << "Norme L2 calculée en fournissant le volume : " << myField1->normL2(myField1_vol) << endl;
 226         for (int i=1; i<=myField1->getNumberOfComponents(); ++i)
 227             cout << "Norme L2 - comp=" << i << " : " << myField1->normL2(i,myField1_vol) << endl;
 228         cout << "Norme L1 calculée en fournissant le volume : " << myField1->normL1(myField1_vol) << endl;
 229         for (int i=1; i<=myField1->getNumberOfComponents(); ++i)
 230             cout << "Norme L1 - comp=" << i << " : " << myField1->normL1(i,myField1_vol) << endl;
 231         delete myField1_vol;
 232
 233         affiche_fieldT(myField1, myField1->getSupport());
 234         cout <<  endl << string(60,'-') << endl;
 235         affiche_fieldT(myField2, myField2->getSupport());
 236         cout << endl << string(60,'-') << endl;
 237
 238         FIELD<double>* myFieldDot = FIELD<double>::scalarProduct(*myField1, *myField2);
 239         affiche_fieldT(myFieldDot, myFieldDot->getSupport());
 240         delete myFieldDot;
 241         cout <<  endl << string(60,'-') << endl ;
 242         affiche_fieldT(&myFieldPlus, myFieldPlus.getSupport());
 243         cout <<  endl << string(60,'-') << endl << endl ;
 244     }
 245
 246
 247     // Verifie plusieurs cas de non compatibilité
 248
 249     // test 1 : Unites non compatibles
 250     const string unite=myField1->getMEDComponentUnit(1);
 251     myField1->setMEDComponentUnit(1,string("UniteBidon"));
 252     ntest++; res=1;
 253     try
 254     {
 255         FIELD<double> myFieldPlus = *myField1 + *myField2;
 256         if(verbose)
 257         {
 258             cout << endl << string(60,'-') << endl;
 259             cout<< "Test " << ntest << " : incompatibilité d'unité : " << endl << endl;
 260         }
 261     }
 262     catch (MEDEXCEPTION & ex)
 263     {
 264         res=0;
 265         if(verbose)
 266             cout << ex.what() << endl;
 267         myField1->setMEDComponentUnit(1,unite);
 268     }
 269     cout << res << endl;
 270
 271     // test 2 : numberOfComponents non compatibles
 272     const int numberOfComponents =myField1->getNumberOfComponents();
 273     myField1->setNumberOfComponents(13);
 274     ntest++; res=1;
 275     try
 276     {
 277         if(verbose)
 278         {
 279             cout << endl << string(60,'-') << endl;
 280             cout<< "Test " << ntest << " : incompatibilité nombre de composantes : " << endl << endl;
 281         }
 282         FIELD<double> myFieldPlus = *myField1 + *myField2;
 283     }
 284     catch (MEDEXCEPTION & ex)
 285     {
 286         res=0;
 287         if(verbose)
 288             cout << endl << ex.what() << endl << endl;
 289         myField1->setNumberOfComponents(numberOfComponents);
 290     }
 291     cout << res << endl;
 292
 293     // test 3 : supports non compatibles
 294     const SUPPORT mySupport2(myMesh,"On_all_node",MED_NODE);
 295     myField1->setSupport(&mySupport2);
 296     ntest++; res=1;
 297     try
 298     {
 299         if(verbose)
 300             cout << endl << string(60,'-') << endl << "Test " << ntest << " : incompatibilité des supports"  << endl << endl;
 301         FIELD<double> myFieldPlus = *myField1 + *myField2;
 302     }
 303     catch (MEDEXCEPTION & ex)
 304     {
 305         res=0;
 306         if(verbose)
 307             cout << ex.what() << endl << endl << endl;
 308         myField1->setSupport(mySupport);
 309     }
 310     cout << res << endl;
 311
 312     // test 4 : champs de taille nulle
 313     myField1->setNumberOfComponents(0);
 314     myField2->setNumberOfComponents(0);
 315     ntest++; res=2;
 316     try
 317     {
 318         if(verbose)
 319             cout<< endl << string(60,'-') << endl << "Test " << ntest << " : incompatibilité taille nulle" << endl << endl;
 320         FIELD<double> myFieldPlus = *myField1 + *myField2;
 321     }
 322     catch (MEDEXCEPTION & ex)
 323     {
 324         --res;
 325         if(verbose)
 326             cout << ex.what() << endl << endl ;
 327     }
 328     try
 329     {
 330         double mynorm2=myField1->norm2();
 331     }
 332     catch (MEDEXCEPTION & ex)
 333     {
 334         --res;
 335         if(verbose)
 336             cout << ex.what() << endl << endl ;
 337         myField1->setNumberOfComponents(numberOfComponents);
 338         myField2->setNumberOfComponents(numberOfComponents);
 339     }
 340     cout << res << endl;
 341
 342     // Apres toutes ces exceptions, des opérations qui marchent!
 343
 344     if(verbose)
 345     {
 346         cout<< endl << string(60,'-') << endl << "Test " << ++ntest << " : Operations arithmétiques" << endl;
 347         cout << endl << " f1           : "; affiche_valeur_field(*myField1);
 348         cout << endl << " f2           : "; affiche_valeur_field(*myField2);
 349         cout  << endl << string(140,'-');
 350     }
 351
 352     // Test du résultats de certaines opérations et affichage si verbose
 353     checkOperation(myFieldPlus, *myField1, *myField2, '+', " f1+f2    ", verbose);
 354     FIELD<double>* myFieldadd = FIELD<double>::add(*myField1, *myField2);
 355     checkOperation( *myFieldadd, *myField1, *myField2, '+', "add(f1,f2)", verbose);
 356     delete myFieldadd;
 357
 358     FIELD<double> myFieldMoins = *myField1 - *myField2;
 359     checkOperation(myFieldMoins, *myField1, *myField2, '-', " f1-f2    ", verbose);
 360     FIELD<double>* myFieldsub = FIELD<double>::sub(*myField1, *myField2);
 361     checkOperation( *myFieldsub, *myField1, *myField2, '-', "sub(f1,f2)", verbose);
 362     delete myFieldsub;
 363     FIELD<double> myFieldNeg = -(*myField1);
 364     checkOperation(myFieldNeg, *myField1, *myField1, 'n', " -f1      ", verbose);
 365
 366     FIELD<double> myFieldFois = *myField1 * *myField2;
 367     checkOperation(myFieldFois, *myField1, *myField2, '*', " f1*f2    ", verbose);
 368     FIELD<double>* myFieldmul = FIELD<double>::mul(*myField1, *myField2);
 369     checkOperation( *myFieldmul, *myField1, *myField2, '*', "mul(f1,f2)", verbose);
 370
 371     FIELD<double> myFieldDiv = *myField1 / *myField2;
 372     checkOperation(myFieldDiv, *myField1, *myField2, '/', " f1/f2    ", verbose);
 373     FIELD<double>* myFielddiv = FIELD<double>::div(*myField1, *myField2);
 374     checkOperation( *myFielddiv, *myField1, *myField2, '/', "div(f1,f2)", verbose);
 375     delete myFielddiv;
 376
 377     FIELD<double> myFieldAsso = (*myField1)+(*myField2)*(*myField2);
 378     checkOperation(myFieldAsso, *myField1, *myField2, 'a', " f1+f2*f2 ", verbose);
 379
 380     myField1->applyLin(4.0,1.0);
 381     checkOperation(*myField1, *myField2, *myField2, 'l', " 4.f1 + 1 ", verbose);
 382     myField1->applyFunc<myfunction1>();
 383     checkOperation( *myField1, *myField2, *myField1, '=', "CB : ->f1)", verbose);
 384
 385     *myField1 += *myField2;
 386     checkOperation(*myField1, *myField2, *myField2, '+', " f1+=f2   ", verbose);
 387
 388     *myField1 -= *myField2;
 389     checkOperation(*myField1, *myField2, *myField2, '=', " f1-=f2   ", verbose);
 390
 391     *myField1 *= *myField2;
 392     checkOperation(*myField1, *myField2, *myField2, '*', " f1*=f2   ", verbose);
 393     *myField1 /= *myField2;
 394     checkOperation(*myField1, *myFieldmul, *myField2, '/', " f1/=f2   ", verbose);
 395     delete myFieldmul;
 396
 397
 398     delete myField1;
 399     delete myField2;
 400     delete mySupport ;
 401     delete myMesh ;
 402     return 0;
 403 }