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