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