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