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