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