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