src/MEDMEMBinTest/test_operation_fieldint.cxx

   1 // Copyright (C) 2007-2012  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 using namespace MEDMEM;
  39 using namespace MED_EN;
  40
  41 int myfunction1(int x);
  42 int myfunction1(int x)
  43 {
  44   return 2*x;
  45 }
  46
  47 int myfunction2(int x);
  48 int myfunction2(int x)
  49 {
  50   return x/2;
  51 }
  52
  53 using namespace std;
  54 static void affiche_field_(FIELD_ * myField, const SUPPORT * mySupport)
  55 {
  56   cout << "Field "<< myField->getName() << " : " <<myField->getDescription() <<  endl ;
  57   int NumberOfComponents = myField->getNumberOfComponents() ;
  58   cout << "- Nombre de composantes : "<< NumberOfComponents << endl ;
  59   cout << "- Nombre de valeurs     : "<< myField->getNumberOfValues() << endl ;
  60   for (int i=1; i<NumberOfComponents+1; i++)
  61     {
  62       cout << "  - composante "<<i<<" :"<<endl ;
  63       cout << "      - nom         : "<<myField->getComponentName(i)<< endl;
  64       cout << "      - description : "<<myField->getComponentDescription(i) << endl;
  65       cout << "      - units       : "<<myField->getMEDComponentUnit(i) << endl;
  66     }
  67   cout << "- iteration :" << endl ;
  68   cout << "    - numero : " << myField->getIterationNumber()<< endl  ;
  69   cout << "    - ordre  : " << myField->getOrderNumber()<< endl  ;
  70   cout << "    - temps  : " << myField->getTime()<< endl  ;
  71
  72   cout << "- Type : " << myField->getValueType()<< endl;
  73
  74   cout << "- Adresse support : " << mySupport << endl;
  75 }
  76
  77 static void affiche_fieldT(FIELD<int> * myField, const SUPPORT * mySupport)
  78 {
  79   affiche_field_((FIELD_ *) myField, mySupport);
  80
  81   cout << "- Valeurs :"<<endl;
  82   int NumberOf = mySupport->getNumberOfElements(MED_ALL_ELEMENTS);
  83   int NumberOfComponents = myField->getNumberOfComponents() ;
  84
  85   for (int i=1; i<NumberOf+1; i++)
  86     {
  87       const int * value = myField->getRow(i) ;
  88       for (int j=0; j<NumberOfComponents; j++)
  89         cout << value[j]<< " ";
  90       cout<<endl;
  91     }
  92   std::cout << std::endl;
  93   std::cout << "Norme euclidienne : " << myField->norm2() << endl;
  94   std::cout << "Norme max         : " << myField->normMax() << endl;
  95   try
  96     {
  97       for (int i=1; i<=myField->getNumberOfComponents(); ++i)
  98         std::cout << "Norme L2 - comp=" << i << " : " << myField->normL2(i) << endl;
  99       std::cout << "Norme L2          : " << myField->normL2() << endl;
 100
 101       for (int i=1; i<=myField->getNumberOfComponents(); ++i)
 102         std::cout << "Norme L1 - comp=" << i << " : " << myField->normL1(i) << endl;
 103       std::cout << "Norme L1          : " << myField->normL1() << endl;
 104     }
 105   catch (MEDEXCEPTION &ex)
 106     {
 107       std::cout << ex.what() << std::endl;
 108     }
 109 }
 110
 111 static void affiche_valeur_field(const char * intitule, const int taille, const FIELD<int>& f)
 112 {
 113   const int * value=f.getValue();
 114   std::cout << endl << intitule;
 115   for(int i=0;i<taille;i++)
 116     std::cout << setw(3) << value[i] << " ";
 117 }
 118
 119 int main (int argc, char ** argv)
 120 {
 121   if (argc != 4)
 122     {
 123       cerr << "Usage : " << argv[0]
 124            << " filename meshname fieldname" << endl << endl;
 125       exit(-1);
 126     }
 127   string filename = argv[1] ;
 128   string meshname = argv[2] ;
 129   string fieldname = argv[3];
 130
 131   MESH * myMesh = new MESH(MED_DRIVER,filename,meshname);
 132   const SUPPORT * mySupport;
 133   FIELD<int> * myField1;
 134   try
 135     {
 136       /* read MESH, SUPPORT and FIELD */
 137       mySupport = myMesh->getSupportOnAll(MED_CELL);
 138       myField1 = new FIELD<int>(mySupport,MED_DRIVER,filename,fieldname) ;
 139     }
 140   catch (MEDEXCEPTION &ex)
 141     {
 142       mySupport = myMesh->getSupportOnAll(MED_NODE);
 143       try
 144         {
 145           myField1 = new FIELD<int>(mySupport,MED_DRIVER,filename,fieldname) ;
 146           myField1->setValueIJ(10,1,-9); // pour tester les normes max avec une valeur negative
 147         }
 148       catch (...)
 149         {
 150           cout << "Field int " << fieldname << " not found !!!" << endl ;
 151           exit (-1) ;
 152         }
 153     }
 154
 155     FIELD<int> * myField2 = new FIELD<int>(* myField1);
 156     //myField1->setNumberOfValues(16); // PROVISOIRE !! BUG
 157     //myField2->setNumberOfValues(16); // PROVISOIRE !! BUG
 158 //      FIELD<int>* myField1_vol=myField1->getSupport()->getMesh()->getVolume(myField1->getSupport());
 159 //      affiche_fieldT(myField1_vol, myField1->getSupport());
 160
 161     affiche_fieldT(myField1, myField1->getSupport());
 162     std::cout <<  endl << string(60,'-') << endl;
 163     affiche_fieldT(myField2, myField2->getSupport());
 164
 165     // Verifie plusieurs cas de non compatibilité
 166
 167     // Unites non compatibles
 168     const string unite=myField1->getMEDComponentUnit(1);
 169     myField1->setMEDComponentUnit(1,string("UniteBidon"));
 170     try
 171     {
 172         std::cout << endl << string(60,'-') << endl;
 173         std::cout<< "Test incompatibilité d'unité :" << endl;
 174         FIELD<int> *myFieldPlus = *myField1 + *myField2;
 175         myFieldPlus->removeReference();
 176     }
 177     catch (MEDEXCEPTION & ex)
 178     {
 179         std::cout << "MEDEXCEPTION : " << ex.what() << endl;
 180         myField1->setMEDComponentUnit(1,unite);
 181     }
 182
 183     // numberOfComponents non compatibles
 184     const int numberOfComponents =myField1->getNumberOfComponents();
 185     myField1->setNumberOfComponents(4);
 186     try
 187     {
 188         std::cout << endl << string(60,'-') << endl;
 189         std::cout<< "Test incompatibilité nombre de composantes :" << endl;
 190         FIELD<int> *myFieldPlus = *myField1 + *myField2;
 191         myFieldPlus->removeReference();
 192     }
 193     catch (MEDEXCEPTION & ex)
 194     {
 195         std::cout << ex.what() << endl;
 196         myField1->setNumberOfComponents(numberOfComponents);
 197     }
 198
 199     // supports non compatibles
 200     const SUPPORT *mySupport2=myMesh->getSupportOnAll(MED_NODE);
 201     myField1->setSupport(mySupport2);
 202     try
 203     {
 204         std::cout << endl << string(60,'-') << endl;
 205         std::cout<< "Test incompatibilité des supports :" << endl;
 206         FIELD<int> *myFieldPlus = *myField1 + *myField2;
 207         myFieldPlus->removeReference();
 208     }
 209     catch (MEDEXCEPTION & ex)
 210     {
 211         std::cout << ex.what() << endl;
 212         myField1->setSupport( myField2->getSupport() );
 213     }
 214
 215     // champs de taille nulle
 216     myField1->setNumberOfComponents(0);
 217     myField2->setNumberOfComponents(0);
 218     try
 219     {
 220         std::cout << endl << string(60,'-') << endl;
 221         std::cout<< "Test incompatibilité taille nulle :" << endl;
 222         FIELD<int> *myFieldPlus = *myField1 + *myField2;
 223         myFieldPlus->removeReference();
 224     }
 225     catch (MEDEXCEPTION & ex)
 226     {
 227         std::cout << ex.what() << endl;
 228     }
 229     try
 230     {
 231         myField1->norm2();
 232     }
 233     catch (MEDEXCEPTION & ex)
 234     {
 235         std::cout << ex.what() << endl;
 236         myField1->setNumberOfComponents(numberOfComponents);
 237         myField2->setNumberOfComponents(numberOfComponents);
 238     }
 239
 240     // Apres toutes ces exceptions, des opérations qui marchent!
 241
 242     FIELD<int> *myFieldPlus = *myField1 + *myField2;
 243     FIELD<int> *myFieldMoins = *myField1 - *myField2;
 244     FIELD<int> *myFieldNeg = -(*myField1);
 245     FIELD<int> *myFieldFois = *myField1 * *myField2;
 246     FIELD<int> *myFieldDiv = *myField1 / *myField2;
 247     FIELD<int> *myFieldAsso = (*myField1)+*((*myField2)*(*myField2));
 248     FIELD<int>* myFieldadd = FIELD<int>::add(*myField1, *myField2);
 249     FIELD<int>* myFieldsub = FIELD<int>::sub(*myField1, *myField2);
 250     FIELD<int>* myFieldmul = FIELD<int>::mul(*myField1, *myField2);
 251     FIELD<int>* myFielddiv = FIELD<int>::div(*myField1, *myField2);
 252     FIELD<int>* myFieldDot = FIELD<int>::scalarProduct(*myField1, *myField2);
 253
 254     std::cout <<  endl << string(60,'-') << endl << "f1+f2 :" << endl << endl;
 255     affiche_fieldT(myFieldPlus, myFieldPlus->getSupport());
 256     std::cout <<  endl << string(60,'-') << endl << "add(f1,f2) :" << endl << endl;
 257     affiche_fieldT(myFieldadd, myFieldadd->getSupport());
 258     std::cout <<  endl << string(60,'-') << endl << "scalarProduct(f1,f2) :" << endl << endl;
 259     affiche_fieldT(myFieldDot, myFieldDot->getSupport());
 260     std::cout <<  endl << string(60,'-') << endl << " - f1 :" << endl << endl;
 261     affiche_fieldT(myFieldNeg, myFieldNeg->getSupport());
 262     int size=myFieldPlus->getNumberOfValues()*myFieldPlus->getNumberOfComponents();
 263
 264     std::cout <<  endl << string(60,'-') << endl << "Tests opérations :" << endl << endl;
 265     affiche_valeur_field("  f1    :", size, *myField1);
 266     affiche_valeur_field("  f2    :", size, *myField2);
 267     std::cout << endl << "        " << string(4*size,'-');
 268
 269     affiche_valeur_field("  +     :", size, *myFieldPlus);
 270     affiche_valeur_field(" add    :", size, *myFieldadd);
 271     affiche_valeur_field("  -     :", size, *myFieldMoins);
 272     affiche_valeur_field(" sub    :", size, *myFieldsub);
 273     affiche_valeur_field("  *     :", size, *myFieldFois);
 274     affiche_valeur_field(" mul    :", size, *myFieldmul);
 275     affiche_valeur_field("  /     :", size, *myFieldDiv);
 276     affiche_valeur_field(" div    :", size, *myFielddiv);
 277     affiche_valeur_field("f1+f2*f1:", size, *myFieldAsso);
 278     affiche_valeur_field("  - f1  :", size, *myFieldNeg);
 279
 280     myFieldPlus->removeReference();
 281     myFieldMoins->removeReference();
 282     myFieldFois->removeReference();
 283     myFieldDiv->removeReference();
 284     myFieldAsso->removeReference();
 285     myFieldNeg->removeReference();
 286
 287     // Test applyLin
 288     std::cout << endl;
 289     myField1->applyLin(1,1);
 290     affiche_valeur_field(" f1+1 :", size, *myField1);
 291     myField1->applyLin(1,-1);
 292     affiche_valeur_field(" -> f1  :", size, *myField1);
 293
 294     // Test applyFunc
 295     std::cout << endl;
 296     myField1->applyFunc<myfunction1>();
 297     affiche_valeur_field(" CB 2f1 :", size, *myField1);
 298     myField1->applyFunc<myfunction2>();
 299     affiche_valeur_field(" -> f1  :", size, *myField1);
 300
 301     // Test operateur +=
 302     std::cout << endl;
 303     *myField1 += *myField2;
 304     affiche_valeur_field(" f1+=f2 :", size, *myField1);
 305
 306     // Test operateur *=
 307     *myField1 *= *myField2;
 308     affiche_valeur_field(" f1*=f2 :", size, *myField1);
 309
 310     // Test operateur /=
 311     *myField1 /= *myField2;
 312     affiche_valeur_field(" f1/=f2 :", size, *myField1);
 313
 314     // Test operateur -=
 315     *myField1 -= *myField2;
 316     affiche_valeur_field(" f1-=f2 :", size, *myField1);
 317
 318     std::cout << endl << endl;
 319
 320
 321     myFieldadd->removeReference();
 322     myFieldsub->removeReference();
 323     myFieldmul->removeReference();
 324     myFielddiv->removeReference();
 325     myFieldDot->removeReference();
 326 //    myField1_vol->removeReference();
 327
 328     myField1->removeReference();
 329     myField2->removeReference();
 330     mySupport->removeReference();
 331     myMesh->removeReference();
 332     return 0;
 333 }