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