src/MEDMEMBinTest/test_operation_fielddouble.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22 #include "MEDMEM_Exception.hxx"
  23 #include "MEDMEM_Mesh.hxx"
  24 #include "MEDMEM_Family.hxx"
  25 #include "MEDMEM_Group.hxx"
  26
  27 #include "MEDMEM_MedMeshDriver.hxx"
  28 #include "MEDMEM_MedFieldDriver.hxx"
  29 #include "MEDMEM_Support.hxx"
  30 #include "MEDMEM_Field.hxx"
  31 #include "MEDMEM_define.hxx"
  32
  33 #include <string>
  34 #include <iostream>
  35 #include <iomanip>
  36 #include <cmath>
  37
  38 double myfunction1(double x)
  39 {
  40   return 0.25*(x-1.0);
  41 }
  42
  43
  44 using namespace std;
  45 using namespace MEDMEM;
  46 using namespace MED_EN;
  47
  48 static void affiche_field_(FIELD_ * myField, const SUPPORT * mySupport)
  49 {
  50   cout << "Field "<< myField->getName() << " : " <<myField->getDescription() <<  endl ;
  51   int NumberOfComponents = myField->getNumberOfComponents() ;
  52   cout << "- Nombre de composantes : "<< NumberOfComponents << endl ;
  53   cout << "- Nombre de valeurs     : "<< myField->getNumberOfValues() << endl ;
  54   for (int i=1; i<NumberOfComponents+1; i++)
  55     {
  56       cout << "  - composante "<<i<<" :"<<endl ;
  57       cout << "      - nom         : "<<myField->getComponentName(i)<< endl;
  58       cout << "      - description : "<<myField->getComponentDescription(i) << endl;
  59       cout << "      - units       : "<<myField->getMEDComponentUnit(i) << endl;
  60     }
  61   cout << "- iteration :" << endl ;
  62   cout << "    - numero : " << myField->getIterationNumber()<< endl  ;
  63   cout << "    - ordre  : " << myField->getOrderNumber()<< endl  ;
  64   cout << "    - temps  : " << myField->getTime()<< endl  ;
  65
  66   cout << "- Type : " << myField->getValueType()<< endl;
  67
  68   cout << "- Adresse support : " << mySupport << endl;
  69 }
  70
  71 static void affiche_fieldT(FIELD<double> * myField, const SUPPORT * mySupport)
  72 {
  73   affiche_field_((FIELD_ *) myField, mySupport);
  74
  75   cout << "- Valeurs :"<<endl;
  76   int NumberOf = mySupport->getNumberOfElements(MED_ALL_ELEMENTS);
  77   int NumberOfComponents = myField->getNumberOfComponents() ;
  78
  79   for (int i=1; i<NumberOf+1; i++)
  80     {
  81       const double * value = myField->getRow(i) ;
  82       for (int j=0; j<NumberOfComponents; j++)
  83         cout << value[j]<< " ";
  84       cout<<endl;
  85     }
  86   cout << endl;
  87   cout << "Norme euclidienne : " << myField->norm2() << endl;
  88   cout << "Norme max         : " << myField->normMax() << endl;
  89   try
  90     {
  91       for (int i=1; i<=myField->getNumberOfComponents(); ++i)
  92         cout << "Norme L2 - comp=" << i << " : " << myField->normL2(i) << endl;
  93       cout << "Norme L2          : " << myField->normL2() << endl;
  94
  95       for (int i=1; i<=myField->getNumberOfComponents(); ++i)
  96         cout << "Norme L1 - comp=" << i << " : " << myField->normL1(i) << endl;
  97       cout << "Norme L1          : " << myField->normL1() << endl;
  98     }
  99   catch (MEDEXCEPTION &ex)
 100     {
 101       cout << ex.what() << endl;
 102     }
 103 }
 104
 105 static void affiche_valeur_field(const FIELD<double>& f)
 106 {
 107   const int tailleMax=12;
 108   const int taille=f.getNumberOfValues()*f.getNumberOfComponents();
 109   const double * value=f.getValue();
 110   if(taille<=tailleMax)
 111     for(int i=0;i<taille;i++)
 112       cout << setw(3) << value[i] << " ";
 113   else
 114     {
 115       for(int i=0; i<tailleMax/2; ++i)
 116         cout << setw(3) << value[i] << " ";
 117       cout << "    ...    ";
 118       for(int i=taille-tailleMax/2 ; i<taille; ++i)
 119         cout << setw(3) << value[i] << " ";
 120     }
 121 }
 122
 123 static void checkOperation(const FIELD<double>& resOp, const FIELD<double>& f1, const FIELD<double>& f2,
 124                            char Op, const char* intitule, int verbose)
 125 {
 126   int res=0;
 127
 128   // get pointers to inside arrays of values
 129   const double * value=resOp.getValue();
 130   const double * value1=f1.getValue();
 131   const double * value2=f2.getValue();
 132   const int size=f1.getNumberOfValues()*f1.getNumberOfComponents(); // size of field1
 133
 134   // check size compatibility
 135   if(f1.getNumberOfValues()*f1.getNumberOfComponents()!=size ||
 136      resOp.getNumberOfValues()*resOp.getNumberOfComponents()!=size)
 137     res=1;
 138
 139   if(!res)
 140     {
 141       switch(Op)
 142         {
 143         case '+':
 144           for(int i=0; i!=size; ++i)
 145             if(value[i]!=value1[i]+value2[i])
 146               res+=1;
 147           break;
 148         case '-':
 149           for(int i=0; i!=size; ++i)
 150             if(value[i]!=value1[i]-value2[i])
 151               res+=1;
 152           break;
 153         case 'n':
 154           for(int i=0; i!=size; ++i)
 155             if(value[i]!=-value1[i])
 156               res+=1;
 157           break;
 158         case '*':
 159           for(int i=0; i!=size; ++i)
 160             if(value[i]!=value1[i]*value2[i])
 161               res+=1;
 162           break;
 163         case '/':
 164           for(int i=0; i!=size; ++i)
 165             if(value2[i]!=0.0)
 166               if(value[i]!=value1[i]/value2[i])
 167                 res+=1;
 168           break;
 169         case '=':
 170           for(int i=0; i!=size; ++i)
 171             if(value[i]!=value2[i])
 172               res+=1;
 173           break;
 174         case 'a':
 175           for(int i=0; i!=size; ++i)
 176             if(value[i]!=value1[i]+value2[i]*value2[i])
 177               res+=1;
 178           break;
 179         }
 180     }
 181
 182   if (verbose)
 183     cout << endl << intitule << "[";
 184   cout << res;
 185   if (verbose)
 186     {
 187       cout << "] : ";
 188       affiche_valeur_field(resOp);
 189     }
 190   else
 191     cout << endl;
 192 }
 193
 194 int main (int argc, char ** argv)
 195 {
 196   /* process the arguments */
 197   int verbose=0;  //  verbose=1 if the verbose mode is selected
 198   int res=0; // unit test result
 199   int ntest=0;  // numÃ©ro du test
 200
 201   if (argc>=2 && !strcmp(argv[1],"-v"))
 202     verbose=1;
 203
 204   if (argc != 4+verbose)
 205     {
 206       cerr << "Usage : " << argv[0]
 207            << "[-v] filename meshname fieldname" << endl << endl
 208            << "-> tests field's operations on the FIELD<double> fieldname" << endl
 209            << "Use optional option -v to select verbose mode" << endl;
 210       exit(-1);
 211     }
 212   string filename  = argv[verbose+1];
 213   string meshname  = argv[verbose+2];// Maintenant plus très utile
 214   string fieldname = argv[verbose+3];
 215
 216   /* read MESH, SUPPORT and FIELDS */
 217   //MESH * myMesh = new MESH(MED_DRIVER,filename,meshname);
 218
 219   MESH * myMesh;
 220   const SUPPORT * mySupport;
 221   FIELD<double> * myField1;
 222
 223   try
 224     {
 225       myField1  = new FIELD<double>(MED_DRIVER,filename,fieldname) ;
 226       mySupport = myField1->getSupport();
 227       myMesh    = new MESH(MED_DRIVER,filename,mySupport->getMeshName());
 228       mySupport->setMesh(myMesh);
 229
 230       FIELD<double> * myField2 = new FIELD<double>(* myField1);
 231       FIELD<double> *myFieldPlus = *myField1 + *myField2;
 232       if(verbose)
 233         {
 234           // affichage des nprmes,des champs f1, f2, scalarProduct(f1,f2) et f1+f2
 235           FIELD<double>* myField1_vol=myField1->getSupport()->getMesh()->getVolume(myField1->getSupport());
 236           cout << "Norme L2 calculee en fournissant le volume : " << myField1->normL2(myField1_vol) << endl;
 237           for (int i=1; i<=myField1->getNumberOfComponents(); ++i)
 238             cout << "Norme L2 - comp=" << i << " : " << myField1->normL2(i,myField1_vol) << endl;
 239           cout << "Norme L1 calculee en fournissant le volume : " << myField1->normL1(myField1_vol) << endl;
 240           for (int i=1; i<=myField1->getNumberOfComponents(); ++i)
 241             cout << "Norme L1 - comp=" << i << " : " << myField1->normL1(i,myField1_vol) << endl;
 242           myField1_vol->removeReference();
 243
 244           affiche_fieldT(myField1, myField1->getSupport());
 245           cout <<  endl << string(60,'-') << endl;
 246           affiche_fieldT(myField2, myField2->getSupport());
 247           cout << endl << string(60,'-') << endl;
 248
 249           FIELD<double>* myFieldDot = FIELD<double>::scalarProduct(*myField1, *myField2);
 250           affiche_fieldT(myFieldDot, myFieldDot->getSupport());
 251           myFieldDot->removeReference();
 252           cout <<  endl << string(60,'-') << endl ;
 253           affiche_fieldT(myFieldPlus, myFieldPlus->getSupport());
 254           cout <<  endl << string(60,'-') << endl << endl ;
 255         }
 256
 257
 258       // Verifie plusieurs cas de non compatibilitÃ©
 259
 260       // test 1 : Unites non compatibles
 261       const string unite=myField1->getMEDComponentUnit(1);
 262       myField1->setMEDComponentUnit(1,string("UniteBidon"));
 263       ntest++; res=1;
 264       try
 265         {
 266           FIELD<double> *myFieldPlus = *myField1 + *myField2;
 267           if(verbose)
 268             {
 269               cout << endl << string(60,'-') << endl;
 270               cout<< "Test " << ntest << " : incompatibilitÃ© d'unitÃ© : " << endl << endl;
 271             }
 272           myFieldPlus->removeReference();
 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           myFieldPlus->removeReference();
 296         }
 297       catch (MEDEXCEPTION & ex)
 298         {
 299           res=0;
 300           if(verbose)
 301             cout << endl << ex.what() << endl << endl;
 302           myField1->setNumberOfComponents(numberOfComponents);
 303         }
 304       cout << res << endl;
 305
 306       // test 3 : supports non compatibles
 307       const SUPPORT *mySupport2=myMesh->getSupportOnAll(MED_NODE);
 308       myField1->setSupport(mySupport2);
 309       ntest++; res=1;
 310       try
 311         {
 312           if(verbose)
 313             cout << endl << string(60,'-') << endl << "Test " << ntest << " : incompatibilitÃ© des supports"  << endl << endl;
 314           FIELD<double> *myFieldPlus = *myField1 + *myField2;
 315           myFieldPlus->removeReference();
 316         }
 317       catch (MEDEXCEPTION & ex)
 318         {
 319           res=0;
 320           if(verbose)
 321             cout << ex.what() << endl << endl << endl;
 322           myField1->setSupport(mySupport);
 323         }
 324       cout << res << endl;
 325
 326       // test 4 : champs de taille nulle
 327       myField1->setNumberOfComponents(0);
 328       myField2->setNumberOfComponents(0);
 329       ntest++; res=2;
 330       try
 331         {
 332           if(verbose)
 333             cout<< endl << string(60,'-') << endl << "Test " << ntest << " : incompatibilitÃ© taille nulle" << endl << endl;
 334           FIELD<double> *myFieldPlus = *myField1 + *myField2;
 335           myFieldPlus->removeReference();
 336         }
 337       catch (MEDEXCEPTION & ex)
 338         {
 339           --res;
 340           if(verbose)
 341             cout << ex.what() << endl << endl ;
 342         }
 343       try
 344         {
 345           myField1->norm2();
 346         }
 347       catch (MEDEXCEPTION & ex)
 348         {
 349           --res;
 350           if(verbose)
 351             cout << ex.what() << endl << endl ;
 352           myField1->setNumberOfComponents(numberOfComponents);
 353           myField2->setNumberOfComponents(numberOfComponents);
 354         }
 355       cout << res << endl;
 356
 357       // Apres toutes ces exceptions, des opÃ©rations qui marchent!
 358
 359       if(verbose)
 360         {
 361           cout<< endl << string(60,'-') << endl << "Test " << ++ntest << " : Operations arithmÃ©tiques" << endl;
 362           cout << endl << " f1           : "; affiche_valeur_field(*myField1);
 363           cout << endl << " f2           : "; affiche_valeur_field(*myField2);
 364           cout  << endl << string(140,'-');
 365         }
 366
 367       // Test du rÃ©sultats de certaines opÃ©rations et affichage si verbose
 368       checkOperation(*myFieldPlus, *myField1, *myField2, '+', " f1+f2    ", verbose);
 369       FIELD<double>* myFieldadd = FIELD<double>::add(*myField1, *myField2);
 370       checkOperation( *myFieldadd, *myField1, *myField2, '+', "add(f1,f2)", verbose);
 371       myFieldadd->removeReference();
 372
 373       FIELD<double> *myFieldMoins = *myField1 - *myField2;
 374       checkOperation(*myFieldMoins, *myField1, *myField2, '-', " f1-f2    ", verbose);
 375       myFieldMoins->removeReference();
 376       FIELD<double>* myFieldsub = FIELD<double>::sub(*myField1, *myField2);
 377       checkOperation( *myFieldsub, *myField1, *myField2, '-', "sub(f1,f2)", verbose);
 378       myFieldsub->removeReference();
 379       FIELD<double> *myFieldNeg = -(*myField1);
 380       checkOperation(*myFieldNeg, *myField1, *myField1, 'n', " -f1      ", verbose);
 381       myFieldNeg->removeReference();
 382       FIELD<double> *myFieldFois = *myField1 * *myField2;
 383       checkOperation(*myFieldFois, *myField1, *myField2, '*', " f1*f2    ", verbose);
 384       myFieldFois->removeReference();
 385       FIELD<double>* myFieldmul = FIELD<double>::mul(*myField1, *myField2);
 386       checkOperation( *myFieldmul, *myField1, *myField2, '*', "mul(f1,f2)", verbose);
 387
 388       FIELD<double> *myFieldDiv = *myField1 / *myField2;
 389       checkOperation(*myFieldDiv, *myField1, *myField2, '/', " f1/f2    ", verbose);
 390       myFieldDiv->removeReference();
 391       FIELD<double>* myFielddiv = FIELD<double>::div(*myField1, *myField2);
 392       checkOperation( *myFielddiv, *myField1, *myField2, '/', "div(f1,f2)", verbose);
 393       myFielddiv->removeReference();
 394
 395       FIELD<double> *myFieldAsso = (*myField1)+*((*myField2)*(*myField2));
 396       checkOperation(*myFieldAsso, *myField1, *myField2, 'a', " f1+f2*f2 ", verbose);
 397       myFieldAsso->removeReference();
 398       myField1->applyLin(4.0,1.0);
 399       checkOperation(*myField1, *myField2, *myField2, 'l', " 4.f1 + 1 ", verbose);
 400       myField1->applyFunc<myfunction1>();
 401       checkOperation( *myField1, *myField2, *myField1, '=', "CB : ->f1)", verbose);
 402
 403       *myField1 += *myField2;
 404       checkOperation(*myField1, *myField2, *myField2, '+', " f1+=f2   ", verbose);
 405
 406       *myField1 -= *myField2;
 407       checkOperation(*myField1, *myField2, *myField2, '=', " f1-=f2   ", verbose);
 408
 409       *myField1 *= *myField2;
 410       checkOperation(*myField1, *myField2, *myField2, '*', " f1*=f2   ", verbose);
 411       *myField1 /= *myField2;
 412       checkOperation(*myField1, *myFieldmul, *myField2, '/', " f1/=f2   ", verbose);
 413       myFieldmul->removeReference();
 414
 415
 416       myFieldPlus->removeReference();
 417       myField1->removeReference();
 418       myField2->removeReference();
 419       //mySupport->removeReference(); ???
 420       myMesh->removeReference();
 421     }
 422   catch ( MEDEXCEPTION & ex)
 423     {
 424       cout << ex.what() << endl;
 425     }
 426
 427   return 0;
 428 }