Back to index

salome-med  6.5.0
test_operation_fielddouble.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
00004 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 #include "MEDMEM_Exception.hxx"
00023 #include "MEDMEM_Mesh.hxx"
00024 #include "MEDMEM_Family.hxx"
00025 #include "MEDMEM_Group.hxx"
00026 
00027 #include "MEDMEM_MedMeshDriver.hxx"
00028 #include "MEDMEM_MedFieldDriver.hxx"
00029 #include "MEDMEM_Support.hxx"
00030 #include "MEDMEM_Field.hxx"
00031 #include "MEDMEM_define.hxx"
00032 
00033 #include <string>
00034 #include <iostream>
00035 #include <iomanip>
00036 #include <cmath>
00037 
00038 double myfunction1(double x)
00039 {
00040   return 0.25*(x-1.0);
00041 }
00042 
00043 
00044 using namespace std;
00045 using namespace MEDMEM;
00046 using namespace MED_EN;
00047 
00048 static void affiche_field_(FIELD_ * myField, const SUPPORT * mySupport)
00049 {
00050   cout << "Field "<< myField->getName() << " : " <<myField->getDescription() <<  endl ;
00051   int NumberOfComponents = myField->getNumberOfComponents() ;
00052   cout << "- Nombre de composantes : "<< NumberOfComponents << endl ;
00053   cout << "- Nombre de valeurs     : "<< myField->getNumberOfValues() << endl ;
00054   for (int i=1; i<NumberOfComponents+1; i++) 
00055     {
00056       cout << "  - composante "<<i<<" :"<<endl ;
00057       cout << "      - nom         : "<<myField->getComponentName(i)<< endl;
00058       cout << "      - description : "<<myField->getComponentDescription(i) << endl;
00059       cout << "      - units       : "<<myField->getMEDComponentUnit(i) << endl;
00060     }
00061   cout << "- iteration :" << endl ;
00062   cout << "    - numero : " << myField->getIterationNumber()<< endl  ;
00063   cout << "    - ordre  : " << myField->getOrderNumber()<< endl  ;
00064   cout << "    - temps  : " << myField->getTime()<< endl  ;
00065 
00066   cout << "- Type : " << myField->getValueType()<< endl;
00067 
00068   cout << "- Adresse support : " << mySupport << endl;
00069 }
00070 
00071 static void affiche_fieldT(FIELD<double> * myField, const SUPPORT * mySupport)
00072 {
00073   affiche_field_((FIELD_ *) myField, mySupport);
00074 
00075   cout << "- Valeurs :"<<endl;
00076   int NumberOf = mySupport->getNumberOfElements(MED_ALL_ELEMENTS);
00077   int NumberOfComponents = myField->getNumberOfComponents() ;
00078 
00079   for (int i=1; i<NumberOf+1; i++) 
00080     {
00081       const double * value = myField->getRow(i) ;
00082       for (int j=0; j<NumberOfComponents; j++)
00083         cout << value[j]<< " ";
00084       cout<<endl;
00085     }
00086   cout << endl;
00087   cout << "Norme euclidienne : " << myField->norm2() << endl;
00088   cout << "Norme max         : " << myField->normMax() << endl;
00089   try
00090     {
00091       for (int i=1; i<=myField->getNumberOfComponents(); ++i)
00092         cout << "Norme L2 - comp=" << i << " : " << myField->normL2(i) << endl;
00093       cout << "Norme L2          : " << myField->normL2() << endl;
00094 
00095       for (int i=1; i<=myField->getNumberOfComponents(); ++i)
00096         cout << "Norme L1 - comp=" << i << " : " << myField->normL1(i) << endl;
00097       cout << "Norme L1          : " << myField->normL1() << endl;
00098     }
00099   catch (MEDEXCEPTION &ex)
00100     {
00101       cout << ex.what() << endl;
00102     }
00103 }
00104 
00105 static void affiche_valeur_field(const FIELD<double>& f)
00106 {
00107   const int tailleMax=12;
00108   const int taille=f.getNumberOfValues()*f.getNumberOfComponents();
00109   const double * value=f.getValue();
00110   if(taille<=tailleMax)
00111     for(int i=0;i<taille;i++)
00112       cout << setw(3) << value[i] << " ";
00113   else
00114     {
00115       for(int i=0; i<tailleMax/2; ++i)
00116         cout << setw(3) << value[i] << " ";
00117       cout << "    ...    ";
00118       for(int i=taille-tailleMax/2 ; i<taille; ++i)
00119         cout << setw(3) << value[i] << " ";
00120     }
00121 }
00122 
00123 static void checkOperation(const FIELD<double>& resOp, const FIELD<double>& f1, const FIELD<double>& f2,
00124                            char Op, const char* intitule, int verbose)
00125 {
00126   int res=0;
00127 
00128   // get pointers to inside arrays of values
00129   const double * value=resOp.getValue();
00130   const double * value1=f1.getValue();
00131   const double * value2=f2.getValue();
00132   const int size=f1.getNumberOfValues()*f1.getNumberOfComponents(); // size of field1
00133 
00134   // check size compatibility
00135   if(f1.getNumberOfValues()*f1.getNumberOfComponents()!=size ||
00136      resOp.getNumberOfValues()*resOp.getNumberOfComponents()!=size)
00137     res=1;
00138 
00139   if(!res)
00140     {
00141       switch(Op)
00142         {
00143         case '+':
00144           for(int i=0; i!=size; ++i)
00145             if(value[i]!=value1[i]+value2[i])
00146               res+=1;
00147           break;
00148         case '-':
00149           for(int i=0; i!=size; ++i)
00150             if(value[i]!=value1[i]-value2[i])
00151               res+=1;
00152           break;
00153         case 'n':
00154           for(int i=0; i!=size; ++i)
00155             if(value[i]!=-value1[i])
00156               res+=1;
00157           break;
00158         case '*':
00159           for(int i=0; i!=size; ++i)
00160             if(value[i]!=value1[i]*value2[i])
00161               res+=1;
00162           break;
00163         case '/':
00164           for(int i=0; i!=size; ++i)
00165             if(value2[i]!=0.0)
00166               if(value[i]!=value1[i]/value2[i])
00167                 res+=1;
00168           break;
00169         case '=':
00170           for(int i=0; i!=size; ++i)
00171             if(value[i]!=value2[i])
00172               res+=1;
00173           break;
00174         case 'a':
00175           for(int i=0; i!=size; ++i)
00176             if(value[i]!=value1[i]+value2[i]*value2[i])
00177               res+=1;
00178           break;
00179         }
00180     }
00181 
00182   if (verbose)
00183     cout << endl << intitule << "[";
00184   cout << res;
00185   if (verbose)
00186     {
00187       cout << "] : ";
00188       affiche_valeur_field(resOp);
00189     }
00190   else
00191     cout << endl;
00192 }
00193 
00194 int main (int argc, char ** argv)
00195 {
00196   /* process the arguments */
00197   int verbose=0;  //  verbose=1 if the verbose mode is selected
00198   int res=0; // unit test result
00199   int ntest=0;  // numéro du test
00200 
00201   if (argc>=2 && !strcmp(argv[1],"-v"))
00202     verbose=1;
00203 
00204   if (argc != 4+verbose)
00205     {
00206       cerr << "Usage : " << argv[0]
00207            << "[-v] filename meshname fieldname" << endl << endl
00208            << "-> tests field's operations on the FIELD<double> fieldname" << endl
00209            << "Use optional option -v to select verbose mode" << endl;
00210       exit(-1);
00211     }
00212   string filename  = argv[verbose+1];
00213   string meshname  = argv[verbose+2];// Maintenant plus très utile
00214   string fieldname = argv[verbose+3];
00215 
00216   /* read MESH, SUPPORT and FIELDS */
00217   //MESH * myMesh = new MESH(MED_DRIVER,filename,meshname);
00218 
00219   MESH * myMesh;
00220   const SUPPORT * mySupport;
00221   FIELD<double> * myField1;
00222 
00223   try
00224     {
00225       myField1  = new FIELD<double>(MED_DRIVER,filename,fieldname) ;
00226       mySupport = myField1->getSupport();
00227       myMesh    = new MESH(MED_DRIVER,filename,mySupport->getMeshName());
00228       mySupport->setMesh(myMesh);
00229 
00230       FIELD<double> * myField2 = new FIELD<double>(* myField1);
00231       FIELD<double> *myFieldPlus = *myField1 + *myField2;
00232       if(verbose)
00233         {
00234           // affichage des nprmes,des champs f1, f2, scalarProduct(f1,f2) et f1+f2
00235           FIELD<double>* myField1_vol=myField1->getSupport()->getMesh()->getVolume(myField1->getSupport());
00236           cout << "Norme L2 calculee en fournissant le volume : " << myField1->normL2(myField1_vol) << endl;
00237           for (int i=1; i<=myField1->getNumberOfComponents(); ++i)
00238             cout << "Norme L2 - comp=" << i << " : " << myField1->normL2(i,myField1_vol) << endl;
00239           cout << "Norme L1 calculee en fournissant le volume : " << myField1->normL1(myField1_vol) << endl;
00240           for (int i=1; i<=myField1->getNumberOfComponents(); ++i)
00241             cout << "Norme L1 - comp=" << i << " : " << myField1->normL1(i,myField1_vol) << endl;
00242           myField1_vol->removeReference();
00243 
00244           affiche_fieldT(myField1, myField1->getSupport());
00245           cout <<  endl << string(60,'-') << endl;
00246           affiche_fieldT(myField2, myField2->getSupport());
00247           cout << endl << string(60,'-') << endl;
00248 
00249           FIELD<double>* myFieldDot = FIELD<double>::scalarProduct(*myField1, *myField2);
00250           affiche_fieldT(myFieldDot, myFieldDot->getSupport());
00251           myFieldDot->removeReference();
00252           cout <<  endl << string(60,'-') << endl ;
00253           affiche_fieldT(myFieldPlus, myFieldPlus->getSupport());
00254           cout <<  endl << string(60,'-') << endl << endl ;
00255         }
00256 
00257 
00258       // Verifie plusieurs cas de non compatibilité
00259 
00260       // test 1 : Unites non compatibles
00261       const string unite=myField1->getMEDComponentUnit(1);
00262       myField1->setMEDComponentUnit(1,string("UniteBidon"));
00263       ntest++; res=1;
00264       try
00265         {
00266           FIELD<double> *myFieldPlus = *myField1 + *myField2;
00267           if(verbose)
00268             {
00269               cout << endl << string(60,'-') << endl;
00270               cout<< "Test " << ntest << " : incompatibilité d'unité : " << endl << endl;
00271             }
00272           myFieldPlus->removeReference();
00273         }
00274       catch (MEDEXCEPTION & ex)
00275         {
00276           res=0;
00277           if(verbose)
00278             cout << ex.what() << endl;
00279           myField1->setMEDComponentUnit(1,unite);
00280         }
00281       cout << res << endl;
00282 
00283       // test 2 : numberOfComponents non compatibles
00284       const int numberOfComponents =myField1->getNumberOfComponents();
00285       myField1->setNumberOfComponents(13);
00286       ntest++; res=1;
00287       try
00288         {
00289           if(verbose)
00290             {
00291               cout << endl << string(60,'-') << endl;
00292               cout<< "Test " << ntest << " : incompatibilité nombre de composantes : " << endl << endl;
00293             }
00294           FIELD<double> *myFieldPlus = *myField1 + *myField2;
00295           myFieldPlus->removeReference();
00296         }
00297       catch (MEDEXCEPTION & ex)
00298         {
00299           res=0;
00300           if(verbose)
00301             cout << endl << ex.what() << endl << endl;
00302           myField1->setNumberOfComponents(numberOfComponents);
00303         }
00304       cout << res << endl;
00305 
00306       // test 3 : supports non compatibles
00307       const SUPPORT *mySupport2=myMesh->getSupportOnAll(MED_NODE);
00308       myField1->setSupport(mySupport2);
00309       ntest++; res=1;
00310       try
00311         {
00312           if(verbose)
00313             cout << endl << string(60,'-') << endl << "Test " << ntest << " : incompatibilité des supports"  << endl << endl;
00314           FIELD<double> *myFieldPlus = *myField1 + *myField2;
00315           myFieldPlus->removeReference();
00316         }
00317       catch (MEDEXCEPTION & ex)
00318         {
00319           res=0;
00320           if(verbose)
00321             cout << ex.what() << endl << endl << endl;
00322           myField1->setSupport(mySupport);
00323         }
00324       cout << res << endl;
00325 
00326       // test 4 : champs de taille nulle
00327       myField1->setNumberOfComponents(0);
00328       myField2->setNumberOfComponents(0);
00329       ntest++; res=2;
00330       try
00331         {
00332           if(verbose)
00333             cout<< endl << string(60,'-') << endl << "Test " << ntest << " : incompatibilité taille nulle" << endl << endl;
00334           FIELD<double> *myFieldPlus = *myField1 + *myField2;
00335           myFieldPlus->removeReference();
00336         }
00337       catch (MEDEXCEPTION & ex)
00338         {
00339           --res;
00340           if(verbose)
00341             cout << ex.what() << endl << endl ;
00342         }
00343       try
00344         {
00345           myField1->norm2();
00346         }
00347       catch (MEDEXCEPTION & ex)
00348         {
00349           --res;
00350           if(verbose)
00351             cout << ex.what() << endl << endl ;
00352           myField1->setNumberOfComponents(numberOfComponents);
00353           myField2->setNumberOfComponents(numberOfComponents);
00354         }
00355       cout << res << endl;
00356 
00357       // Apres toutes ces exceptions, des opérations qui marchent!
00358 
00359       if(verbose)
00360         {
00361           cout<< endl << string(60,'-') << endl << "Test " << ++ntest << " : Operations arithmétiques" << endl;
00362           cout << endl << " f1           : "; affiche_valeur_field(*myField1);
00363           cout << endl << " f2           : "; affiche_valeur_field(*myField2);
00364           cout  << endl << string(140,'-');
00365         }
00366 
00367       // Test du résultats de certaines opérations et affichage si verbose
00368       checkOperation(*myFieldPlus, *myField1, *myField2, '+', " f1+f2    ", verbose);
00369       FIELD<double>* myFieldadd = FIELD<double>::add(*myField1, *myField2);
00370       checkOperation( *myFieldadd, *myField1, *myField2, '+', "add(f1,f2)", verbose);
00371       myFieldadd->removeReference();
00372 
00373       FIELD<double> *myFieldMoins = *myField1 - *myField2;
00374       checkOperation(*myFieldMoins, *myField1, *myField2, '-', " f1-f2    ", verbose);
00375       myFieldMoins->removeReference();
00376       FIELD<double>* myFieldsub = FIELD<double>::sub(*myField1, *myField2);
00377       checkOperation( *myFieldsub, *myField1, *myField2, '-', "sub(f1,f2)", verbose);
00378       myFieldsub->removeReference();
00379       FIELD<double> *myFieldNeg = -(*myField1);
00380       checkOperation(*myFieldNeg, *myField1, *myField1, 'n', " -f1      ", verbose);
00381       myFieldNeg->removeReference();
00382       FIELD<double> *myFieldFois = *myField1 * *myField2;
00383       checkOperation(*myFieldFois, *myField1, *myField2, '*', " f1*f2    ", verbose);
00384       myFieldFois->removeReference();
00385       FIELD<double>* myFieldmul = FIELD<double>::mul(*myField1, *myField2);
00386       checkOperation( *myFieldmul, *myField1, *myField2, '*', "mul(f1,f2)", verbose);
00387 
00388       FIELD<double> *myFieldDiv = *myField1 / *myField2;
00389       checkOperation(*myFieldDiv, *myField1, *myField2, '/', " f1/f2    ", verbose);
00390       myFieldDiv->removeReference();
00391       FIELD<double>* myFielddiv = FIELD<double>::div(*myField1, *myField2);
00392       checkOperation( *myFielddiv, *myField1, *myField2, '/', "div(f1,f2)", verbose);
00393       myFielddiv->removeReference();
00394 
00395       FIELD<double> *myFieldAsso = (*myField1)+*((*myField2)*(*myField2));
00396       checkOperation(*myFieldAsso, *myField1, *myField2, 'a', " f1+f2*f2 ", verbose);
00397       myFieldAsso->removeReference();
00398       myField1->applyLin(4.0,1.0);
00399       checkOperation(*myField1, *myField2, *myField2, 'l', " 4.f1 + 1 ", verbose);
00400       myField1->applyFunc<myfunction1>();
00401       checkOperation( *myField1, *myField2, *myField1, '=', "CB : ->f1)", verbose);
00402 
00403       *myField1 += *myField2;
00404       checkOperation(*myField1, *myField2, *myField2, '+', " f1+=f2   ", verbose);
00405 
00406       *myField1 -= *myField2;
00407       checkOperation(*myField1, *myField2, *myField2, '=', " f1-=f2   ", verbose);
00408 
00409       *myField1 *= *myField2;
00410       checkOperation(*myField1, *myField2, *myField2, '*', " f1*=f2   ", verbose);
00411       *myField1 /= *myField2;
00412       checkOperation(*myField1, *myFieldmul, *myField2, '/', " f1/=f2   ", verbose);
00413       myFieldmul->removeReference();
00414 
00415 
00416       myFieldPlus->removeReference();
00417       myField1->removeReference();
00418       myField2->removeReference();
00419       //mySupport->removeReference(); ???
00420       myMesh->removeReference();
00421     }
00422   catch ( MEDEXCEPTION & ex) 
00423     {
00424       cout << ex.what() << endl;
00425     }
00426 
00427   return 0;
00428 }