1 // Programme de test des operations sur les champs
8 #include "MEDMEM_Exception.hxx"
9 #include "MEDMEM_Mesh.hxx"
10 #include "MEDMEM_Family.hxx"
11 #include "MEDMEM_Group.hxx"
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"
20 #include "LocalTraceCollector.hxx"
21 #endif /* ifdef _DEBUG_*/
23 double myfunction1(double x)
30 using namespace MEDMEM;
31 using namespace MED_EN;
33 void affiche_field_(FIELD_ * myField, const SUPPORT * mySupport)
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;
45 cout << "- iteration :" << endl ;
46 cout << " - numero : " << myField->getIterationNumber()<< endl ;
47 cout << " - ordre : " << myField->getOrderNumber()<< endl ;
48 cout << " - temps : " << myField->getTime()<< endl ;
50 cout << "- Type : " << myField->getValueType()<< endl;
52 cout << "- Adresse support : " << mySupport << endl;
55 void affiche_fieldT(FIELD<double> * myField, const SUPPORT * mySupport)
57 affiche_field_((FIELD_ *) myField, mySupport);
59 cout << "- Valeurs :"<<endl;
60 int NumberOf = mySupport->getNumberOfElements(MED_ALL_ELEMENTS);
61 int NumberOfComponents = myField->getNumberOfComponents() ;
63 for (int i=1; i<NumberOf+1; i++) {
64 const double * value = myField->getValueI(MED_FULL_INTERLACE,i) ;
65 for (int j=0; j<NumberOfComponents; j++)
66 cout << value[j]<< " ";
70 cout << "Norme euclidienne : " << myField->norm2() << endl;
71 cout << "Norme max : " << myField->normMax() << endl;
74 for (int i=1; i<=myField->getNumberOfComponents(); ++i)
75 cout << "Norme L2 - comp=" << i << " : " << myField->normL2(i) << endl;
76 cout << "Norme L2 : " << myField->normL2() << endl;
78 for (int i=1; i<=myField->getNumberOfComponents(); ++i)
79 cout << "Norme L1 - comp=" << i << " : " << myField->normL1(i) << endl;
80 cout << "Norme L1 : " << myField->normL1() << endl;
82 catch (MEDEXCEPTION &ex)
84 cout << ex.what() << endl;
88 void affiche_valeur_field(const FIELD<double>& f)
90 const int tailleMax=12;
91 const int taille=f.getNumberOfValues()*f.getNumberOfComponents();
92 const double * value=f.getValue(f.getvalue()->getMode());
94 for(int i=0;i<taille;i++)
95 cout << setw(3) << value[i] << " ";
98 for(int i=0; i<tailleMax/2; ++i)
99 cout << setw(3) << value[i] << " ";
101 for(int i=taille-tailleMax/2 ; i<taille; ++i)
102 cout << setw(3) << value[i] << " ";
106 void checkOperation(const FIELD<double>& resOp, const FIELD<double>& f1, const FIELD<double>& f2,
107 char Op, const char* intitule, int verbose)
111 // get pointers to inside arrays of values
112 medModeSwitch mode=resOp.getvalue()->getMode();
113 const double * value=resOp.getValue(mode);
114 const double * value1=f1.getValue(mode);
115 const double * value2=f2.getValue(mode);
116 const int size=f1.getNumberOfValues()*f1.getNumberOfComponents(); // size of field1
118 // check size compatibility
119 if(f1.getNumberOfValues()*f1.getNumberOfComponents()!=size ||
120 resOp.getNumberOfValues()*resOp.getNumberOfComponents()!=size)
128 for(int i=0; i!=size; ++i)
129 if(value[i]!=value1[i]+value2[i])
133 for(int i=0; i!=size; ++i)
134 if(value[i]!=value1[i]-value2[i])
138 for(int i=0; i!=size; ++i)
139 if(value[i]!=-value1[i])
143 for(int i=0; i!=size; ++i)
144 if(value[i]!=value1[i]*value2[i])
148 for(int i=0; i!=size; ++i)
150 if(value[i]!=value1[i]/value2[i])
154 for(int i=0; i!=size; ++i)
155 if(value[i]!=value2[i])
159 for(int i=0; i!=size; ++i)
160 if(value[i]!=value1[i]+value2[i]*value2[i])
168 cout << endl << intitule << "[";
173 affiche_valeur_field(resOp);
179 int main (int argc, char ** argv)
181 /* process the arguments */
182 int verbose=0; // verbose=1 if the verbose mode is selected
183 int res=0; // unit test result
184 int ntest=0; // numéro du test
186 if (argc>=2 && !strcmp(argv[1],"-v"))
189 if (argc != 4+verbose)
191 cerr << "Usage : " << argv[0]
192 << "[-v] filename meshname fieldname" << endl << endl
193 << "-> tests field's operations on the FIELD<double> fieldname" << endl
194 << "Use optional option -v to select verbose mode" << endl;
199 LocalTraceCollector::instance();
200 #endif /* ifdef _DEBUG_*/
202 string filename = argv[verbose+1];
203 string meshname = argv[verbose+2];
204 string fieldname = argv[verbose+3];
206 /* read MESH, SUPPORT and FIELDS */
207 MESH * myMesh = new MESH(MED_DRIVER,filename,meshname);
209 FIELD<double> * myField1;
212 mySupport = new SUPPORT(myMesh,"Support on all Cells",MED_CELL);
213 myField1 = new FIELD<double>(mySupport,MED_DRIVER,filename,fieldname) ;
214 myField1->setValueType(MED_REEL64);
216 catch (MEDEXCEPTION &ex)
218 // field wasn't found on cells, try on nodes
220 mySupport = new SUPPORT(myMesh,"On_all_node",MED_NODE);
223 myField1 = new FIELD<double>(mySupport,MED_DRIVER,filename,fieldname) ;
227 cout << "Field double " << fieldname << " not found !!!" << endl ;
231 FIELD<double> * myField2 = new FIELD<double>(* myField1);
232 FIELD<double> myFieldPlus = *myField1 + *myField2;
235 // affichage des nprmes,des champs f1, f2, scalarProduct(f1,f2) et f1+f2
236 FIELD<double>* myField1_vol=myField1->getSupport()->getMesh()->getVolume(myField1->getSupport());
237 cout << "Norme L2 calculée en fournissant le volume : " << myField1->normL2(myField1_vol) << endl;
238 for (int i=1; i<=myField1->getNumberOfComponents(); ++i)
239 cout << "Norme L2 - comp=" << i << " : " << myField1->normL2(i,myField1_vol) << endl;
240 cout << "Norme L1 calculée en fournissant le volume : " << myField1->normL1(myField1_vol) << endl;
241 for (int i=1; i<=myField1->getNumberOfComponents(); ++i)
242 cout << "Norme L1 - comp=" << i << " : " << myField1->normL1(i,myField1_vol) << endl;
245 affiche_fieldT(myField1, myField1->getSupport());
246 cout << endl << string(60,'-') << endl;
247 affiche_fieldT(myField2, myField2->getSupport());
248 cout << endl << string(60,'-') << endl;
250 FIELD<double>* myFieldDot = FIELD<double>::scalarProduct(*myField1, *myField2);
251 affiche_fieldT(myFieldDot, myFieldDot->getSupport());
253 cout << endl << string(60,'-') << endl ;
254 affiche_fieldT(&myFieldPlus, myFieldPlus.getSupport());
255 cout << endl << string(60,'-') << endl << endl ;
259 // Verifie plusieurs cas de non compatibilité
261 // test 1 : Unites non compatibles
262 const string unite=myField1->getMEDComponentUnit(1);
263 myField1->setMEDComponentUnit(1,string("UniteBidon"));
267 FIELD<double> myFieldPlus = *myField1 + *myField2;
270 cout << endl << string(60,'-') << endl;
271 cout<< "Test " << ntest << " : incompatibilité d'unité : " << endl << endl;
274 catch (MEDEXCEPTION & ex)
278 cout << ex.what() << endl;
279 myField1->setMEDComponentUnit(1,unite);
283 // test 2 : numberOfComponents non compatibles
284 const int numberOfComponents =myField1->getNumberOfComponents();
285 myField1->setNumberOfComponents(13);
291 cout << endl << string(60,'-') << endl;
292 cout<< "Test " << ntest << " : incompatibilité nombre de composantes : " << endl << endl;
294 FIELD<double> myFieldPlus = *myField1 + *myField2;
296 catch (MEDEXCEPTION & ex)
300 cout << endl << ex.what() << endl << endl;
301 myField1->setNumberOfComponents(numberOfComponents);
305 // test 3 : supports non compatibles
306 const SUPPORT mySupport2(myMesh,"On_all_node",MED_NODE);
307 myField1->setSupport(&mySupport2);
312 cout << endl << string(60,'-') << endl << "Test " << ntest << " : incompatibilité des supports" << endl << endl;
313 FIELD<double> myFieldPlus = *myField1 + *myField2;
315 catch (MEDEXCEPTION & ex)
319 cout << ex.what() << endl << endl << endl;
320 myField1->setSupport(mySupport);
324 // test 4 : champs de taille nulle
325 myField1->setNumberOfComponents(0);
326 myField2->setNumberOfComponents(0);
331 cout<< endl << string(60,'-') << endl << "Test " << ntest << " : incompatibilité taille nulle" << endl << endl;
332 FIELD<double> myFieldPlus = *myField1 + *myField2;
334 catch (MEDEXCEPTION & ex)
338 cout << ex.what() << endl << endl ;
342 double mynorm2=myField1->norm2();
344 catch (MEDEXCEPTION & ex)
348 cout << ex.what() << endl << endl ;
349 myField1->setNumberOfComponents(numberOfComponents);
350 myField2->setNumberOfComponents(numberOfComponents);
354 // Apres toutes ces exceptions, des opérations qui marchent!
358 cout<< endl << string(60,'-') << endl << "Test " << ++ntest << " : Operations arithmétiques" << endl;
359 cout << endl << " f1 : "; affiche_valeur_field(*myField1);
360 cout << endl << " f2 : "; affiche_valeur_field(*myField2);
361 cout << endl << string(140,'-');
364 // Test du résultats de certaines opérations et affichage si verbose
365 checkOperation(myFieldPlus, *myField1, *myField2, '+', " f1+f2 ", verbose);
366 FIELD<double>* myFieldadd = FIELD<double>::add(*myField1, *myField2);
367 checkOperation( *myFieldadd, *myField1, *myField2, '+', "add(f1,f2)", verbose);
370 FIELD<double> myFieldMoins = *myField1 - *myField2;
371 checkOperation(myFieldMoins, *myField1, *myField2, '-', " f1-f2 ", verbose);
372 FIELD<double>* myFieldsub = FIELD<double>::sub(*myField1, *myField2);
373 checkOperation( *myFieldsub, *myField1, *myField2, '-', "sub(f1,f2)", verbose);
375 FIELD<double> myFieldNeg = -(*myField1);
376 checkOperation(myFieldNeg, *myField1, *myField1, 'n', " -f1 ", verbose);
378 FIELD<double> myFieldFois = *myField1 * *myField2;
379 checkOperation(myFieldFois, *myField1, *myField2, '*', " f1*f2 ", verbose);
380 FIELD<double>* myFieldmul = FIELD<double>::mul(*myField1, *myField2);
381 checkOperation( *myFieldmul, *myField1, *myField2, '*', "mul(f1,f2)", verbose);
383 FIELD<double> myFieldDiv = *myField1 / *myField2;
384 checkOperation(myFieldDiv, *myField1, *myField2, '/', " f1/f2 ", verbose);
385 FIELD<double>* myFielddiv = FIELD<double>::div(*myField1, *myField2);
386 checkOperation( *myFielddiv, *myField1, *myField2, '/', "div(f1,f2)", verbose);
389 FIELD<double> myFieldAsso = (*myField1)+(*myField2)*(*myField2);
390 checkOperation(myFieldAsso, *myField1, *myField2, 'a', " f1+f2*f2 ", verbose);
392 myField1->applyLin(4.0,1.0);
393 checkOperation(*myField1, *myField2, *myField2, 'l', " 4.f1 + 1 ", verbose);
394 myField1->applyFunc<myfunction1>();
395 checkOperation( *myField1, *myField2, *myField1, '=', "CB : ->f1)", verbose);
397 *myField1 += *myField2;
398 checkOperation(*myField1, *myField2, *myField2, '+', " f1+=f2 ", verbose);
400 *myField1 -= *myField2;
401 checkOperation(*myField1, *myField2, *myField2, '=', " f1-=f2 ", verbose);
403 *myField1 *= *myField2;
404 checkOperation(*myField1, *myField2, *myField2, '*', " f1*=f2 ", verbose);
405 *myField1 /= *myField2;
406 checkOperation(*myField1, *myFieldmul, *myField2, '/', " f1/=f2 ", verbose);