src/MEDMEMCppTest/MEDMEMTest_Field_fault.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "MEDMEMTest.hxx"
  21
  22 #include "MEDMEM_FieldConvert.hxx"
  23 #include "MEDMEM_Field.hxx"
  24 #include "MEDMEM_Mesh.hxx"
  25 #include "MEDMEM_Group.hxx"
  26 #include "MEDMEM_Support.hxx"
  27 #include "MEDMEM_VtkMeshDriver.hxx"
  28 #include "MEDMEM_MedMeshDriver.hxx"
  29
  30 #include <cppunit/TestAssert.h>
  31
  32 #include <sstream>
  33 #include <cmath>
  34
  35 // use this define to enable lines, execution of which leads to Segmentation Fault
  36 //#define ENABLE_FAULTS
  37
  38 // use this define to enable CPPUNIT asserts and fails, showing bugs
  39 //#define ENABLE_FORCED_FAILURES
  40
  41 using namespace std;
  42 using namespace MEDMEM;
  43
  44 // #14,15: MEDMEMTest_Field.cxx
  45 // Check methods from MEDMEM_Field.hxx, MEDMEM_FieldConvert.hxx
  46
  47 /*!
  48  *  Check methods (48), defined in MEDMEM_Field.hxx:
  49  *  class FIELD_
  50 {
  51  *   (+)     FIELD_();
  52  *   (+)     FIELD_(const SUPPORT * Support, const int NumberOfComponents);
  53  *   (+)     FIELD_(const FIELD_ &m);
  54  *   (+)     virtual ~FIELD_();
  55  *   (+)     FIELD_& operator=(const FIELD_ &m);
  56  *
  57  *   (-)     virtual  void     rmDriver(int index=0);
  58  *   (-)     virtual   int     addDriver(driverTypes driverType,
  59  *                             const string & fileName="Default File Name.med",
  60  *                              const string & driverFieldName="Default Field Nam",
  61  *                              MED_EN::med_mode_acces access=MED_EN::MED_REMP);
  62  *   (-)     virtual  int      addDriver(GENDRIVER & driver);
  63  *
  64  *   (-)     virtual  void     read (const GENDRIVER &);
  65  *   (-)     virtual  void     read(int index=0);
  66  *   (-)     virtual  void     openAppend(void);
  67  *   (-)     virtual  void     write(const GENDRIVER &);
  68  *   (-)     virtual  void     write(int index=0, const string & driverName="");
  69  *   (-)     virtual  void     writeAppend(const GENDRIVER &);
  70  *   (-)     virtual  void     writeAppend(int index=0, const string & driverName="");
  71  *
  72  *   (+)     inline void     setName(const string Name);
  73  *   (+)     inline string   getName() const;
  74  *   (+)     inline void     setDescription(const string Description);
  75  *   (+)     inline string   getDescription() const;
  76  *   (+)     inline const SUPPORT * getSupport() const;
  77  *   (+)     inline void     setSupport(const SUPPORT * support);
  78  *   (+)     inline void     setNumberOfComponents(const int NumberOfComponents);
  79  *   (+)     inline int      getNumberOfComponents() const;
  80  *   (+)     inline void     setNumberOfValues(const int NumberOfValues);
  81  *   (+)     inline int      getNumberOfValues() const;
  82  *   (+)     inline void     setComponentsNames(const string * ComponentsNames);
  83  *   (+)     inline void     setComponentName(int i, const string ComponentName);
  84  *   (+)     inline const string * getComponentsNames() const;
  85  *   (+)     inline string   getComponentName(int i) const;
  86  *   (+)     inline void     setComponentsDescriptions(const string * ComponentsDescriptions);
  87  *   (+)     inline void     setComponentDescription(int i, const string ComponentDescription);
  88  *   (+)     inline const string * getComponentsDescriptions() const;
  89  *   (+)     inline string   getComponentDescription(int i) const;
  90  *   (+)     inline void     setComponentsUnits(const UNIT * ComponentsUnits);
  91  *   (+)     inline const UNIT *   getComponentsUnits() const;
  92  *   (+)     inline const UNIT *   getComponentUnit(int i) const;
  93  *   (+)     inline void     setMEDComponentsUnits(const string * MEDComponentsUnits);
  94  *   (+)     inline void     setMEDComponentUnit(int i, const string MEDComponentUnit);
  95  *   (+)     inline const string * getMEDComponentsUnits() const;
  96  *   (+)     inline string   getMEDComponentUnit(int i) const;
  97  *
  98  *   (+)     inline void     setIterationNumber(int IterationNumber);
  99  *   (+)     inline int      getIterationNumber() const;
 100  *   (+)     inline void     setTime(double Time);
 101  *   (+)     inline double   getTime() const;
 102  *   (+)     inline void     setOrderNumber(int OrderNumber);
 103  *   (+)     inline int      getOrderNumber() const;
 104  *
 105  *   (+)     inline MED_EN::med_type_champ getValueType () const;
 106  *   (+)     inline MED_EN::medModeSwitch  getInterlacingType() const;
 107  *   (-)     virtual inline bool getGaussPresence() const throw (MEDEXCEPTION);
 108  *
 109 }
 110  *
 111  *  template <class T, class INTERLACING_TAG> class FIELD : public FIELD_
 112 {
 113  *   (+)     FIELD();
 114  *   (+)     FIELD(const FIELD &m);
 115  *   (+)     FIELD(const SUPPORT * Support, const int NumberOfComponents) throw (MEDEXCEPTION);
 116  *   (+)     FIELD(driverTypes driverType,
 117  *                 const string & fileName, const string & fieldDriverName,
 118  *                 const int iterationNumber=-1, const int orderNumber=-1) throw (MEDEXCEPTION);
 119  *   (+) FIELD(const SUPPORT * Support, driverTypes driverType,
 120  *                 const string & fileName="", const string & fieldName="",
 121  *                 const int iterationNumber = -1, const int orderNumber = -1) throw (MEDEXCEPTION);
 122  *   (+)     ~FIELD();
 123  *   (+)     FIELD & operator=(const FIELD &m);
 124  *
 125  *   (+) const FIELD operator+(const FIELD& m) const;
 126  *   (+) const FIELD operator-(const FIELD& m) const;
 127  *   (+) const FIELD operator*(const FIELD& m) const;
 128  *   (+) const FIELD operator/(const FIELD& m) const;
 129  *   (+) const FIELD operator-() const;
 130  *   (+) FIELD& operator+=(const FIELD& m);
 131  *   (+) FIELD& operator-=(const FIELD& m);
 132  *   (+) FIELD& operator*=(const FIELD& m);
 133  *   (+) FIELD& operator/=(const FIELD& m);
 134  *
 135  *   (+) static FIELD* add(const FIELD& m, const FIELD& n);
 136  *   (+) static FIELD* addDeep(const FIELD& m, const FIELD& n);
 137  *   (+) static FIELD* sub(const FIELD& m, const FIELD& n);
 138  *   (+) static FIELD* subDeep(const FIELD& m, const FIELD& n);
 139  *   (+) static FIELD* mul(const FIELD& m, const FIELD& n);
 140  *   (+) static FIELD* mulDeep(const FIELD& m, const FIELD& n);
 141  *   (+) static FIELD* div(const FIELD& m, const FIELD& n);
 142  *   (+) static FIELD* divDeep(const FIELD& m, const FIELD& n);
 143  *
 144  *   (+)     double normMax() const throw (MEDEXCEPTION);
 145  *   (+)     double norm2() const throw (MEDEXCEPTION);
 146  *
 147  *   (+)     void   applyLin(T a, T b);
 148  *   (+)     template <T T_function(T)> void applyFunc();
 149  *   (+)     void applyPow(T scalar);
 150  *
 151  *   (+)     static FIELD* scalarProduct(const FIELD& m, const FIELD& n, bool deepCheck=false);
 152  *
 153  *   (+)     double normL2(int component, const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 154  *   (+)     double normL2(const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 155  *   (+)     double normL1(int component, const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 156  *   (+)     double normL1(const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 157  *
 158  *   (+)     FIELD* extract(const SUPPORT *subSupport) const throw (MEDEXCEPTION);
 159  *
 160  *   (EMPTY COMMENT, EMPTY IMPLEMENTATION!!!) void init ();
 161  *
 162  *   (+)     void rmDriver(int index=0);
 163  *   (+)     int  addDriver(driverTypes driverType,
 164  *                          const string & fileName="Default File Name.med",
 165  *                          const string & driverFieldName="Default Field Name",
 166  *                          MED_EN::med_mode_acces access=MED_EN::MED_REMP);
 167  *   (+)     int  addDriver(GENDRIVER & driver);
 168  *
 169  *   (+)     void allocValue(const int NumberOfComponents);
 170  *   (+)     void allocValue(const int NumberOfComponents, const int LengthValue);
 171  *   (+)     void deallocValue();
 172  *
 173  *   (+)     inline void read(int index=0);
 174  *   (+)     inline void read(const GENDRIVER & genDriver);
 175  *   (+)     inline void write(int index=0, const string & driverName = "");
 176  *   (+)     inline void write(const GENDRIVER &);
 177  *   (+)     inline void writeAppend(int index=0, const string & driverName = "");
 178  *   (+) inline void writeAppend(const GENDRIVER &);
 179  *
 180  *   (+)     inline MEDMEM_Array_  * getArray()        const throw (MEDEXCEPTION);
 181  *   (+)     inline ArrayGauss     * getArrayGauss()   const throw (MEDEXCEPTION);
 182  *   (+)     inline ArrayNoGauss   * getArrayNoGauss() const throw (MEDEXCEPTION);
 183  *   (+)     inline bool             getGaussPresence() const throw (MEDEXCEPTION);
 184  *
 185  *   (+)     inline int      getValueLength() const throw (MEDEXCEPTION);
 186  *   (+)     inline const T* getValue()       const throw (MEDEXCEPTION);
 187  *   (+)     inline const T* getRow(int i)    const throw (MEDEXCEPTION);
 188  *   (+)     inline const T* getColumn(int j) const throw (MEDEXCEPTION);
 189  *   (+)     inline T        getValueIJ(int i,int j) const throw (MEDEXCEPTION);
 190  *   (+)     inline T        getValueIJK(int i,int j,int k) const throw (MEDEXCEPTION);
 191  *   (+)     bool            getValueOnElement(int eltIdInSup,T* retValues) const throw (MEDEXCEPTION);
 192  *
 193  *   (+)     const int getNumberOfGeometricTypes() const throw (MEDEXCEPTION);
 194  *
 195  *   (+)     const GAUSS_LOCALIZATION<INTERLACING_TAG> & getGaussLocalization
 196  *                          (MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
 197  *   (+)     const GAUSS_LOCALIZATION<INTERLACING_TAG> * getGaussLocalizationPtr
 198  *                          (MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
 199  *   (+)     void setGaussLocalization(MED_EN::medGeometryElement geomElement,
 200  *                                     const GAUSS_LOCALIZATION<INTERLACING_TAG> & gaussloc);
 201  *   (+)     const int * getNumberOfGaussPoints() const throw (MEDEXCEPTION);
 202  *   (+)     const int   getNumberOfGaussPoints
 203  *                          (MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
 204  *   (+)     const int   getNbGaussI(int i) const throw (MEDEXCEPTION);
 205  *
 206  *   (+)     const int * getNumberOfElements() const throw (MEDEXCEPTION);
 207  *   (+)     const MED_EN::medGeometryElement * getGeometricTypes() const throw (MEDEXCEPTION);
 208  *   (+)     bool        isOnAllElements() const throw (MEDEXCEPTION);
 209  *
 210  *   (+)     inline void setArray(MEDMEM_Array_ *value) throw (MEDEXCEPTION);
 211  *   (+)     inline void setValue(T* value) throw (MEDEXCEPTION);
 212  *   (+)     inline void setRow(int i, T* value) throw (MEDEXCEPTION);
 213  *   (+)     inline void setColumn(int i, T* value) throw (MEDEXCEPTION);
 214  *   (+)     inline void setValueIJ(int i, int j, T value) throw (MEDEXCEPTION);
 215  *
 216  *   (NOT IMPLEMENTED!!!) void getVolume() const throw (MEDEXCEPTION);
 217  *   (NOT IMPLEMENTED!!!) void getArea() const throw (MEDEXCEPTION);
 218  *   (NOT IMPLEMENTED!!!) void getLength() const throw (MEDEXCEPTION);
 219  *   (NOT IMPLEMENTED!!!) void getNormal() const throw (MEDEXCEPTION);
 220  *   (NOT IMPLEMENTED!!!) void getBarycenter() const throw (MEDEXCEPTION);
 221  *
 222  *   (+)     void fillFromAnalytic(myFuncType f) throw (MEDEXCEPTION);
 223  *
 224 }
 225  *
 226  *  Use code of test_operation_fieldint.cxx
 227  *              test_operation_fielddouble.cxx
 228  *              test_copie_field_.cxx
 229  *              test_copie_fieldT.cxx
 230  */
 231  static void compareField_(const FIELD_ * theField_1, const FIELD_ * theField_2, bool isFIELD, bool isValue)
 232 {
 233   // name, description, support
 234   CPPUNIT_ASSERT_EQUAL(theField_1->getName(), theField_2->getName());
 235   CPPUNIT_ASSERT_EQUAL(theField_1->getDescription(), theField_2->getDescription());
 236   CPPUNIT_ASSERT_EQUAL(theField_1->getSupport(), theField_2->getSupport());
 237
 238   // components information
 239   int aNbComps = theField_1->getNumberOfComponents();
 240   CPPUNIT_ASSERT_EQUAL(aNbComps, theField_2->getNumberOfComponents());
 241
 242   for (int i = 1; i <= aNbComps; i++)
 243     {
 244       CPPUNIT_ASSERT_EQUAL(theField_1->getComponentName(i), theField_2->getComponentName(i));
 245       CPPUNIT_ASSERT_EQUAL(theField_1->getComponentDescription(i), theField_2->getComponentDescription(i));
 246       CPPUNIT_ASSERT_EQUAL(theField_1->getMEDComponentUnit(i), theField_2->getMEDComponentUnit(i));
 247     }
 248
 249   // iteration information
 250   CPPUNIT_ASSERT_EQUAL(theField_1->getIterationNumber(), theField_2->getIterationNumber());
 251   CPPUNIT_ASSERT_EQUAL(theField_1->getOrderNumber(), theField_2->getOrderNumber());
 252   CPPUNIT_ASSERT_DOUBLES_EQUAL(theField_1->getTime(), theField_2->getTime(), 0.0000001);
 253
 254   // Value
 255   int nbOfValues = theField_1->getNumberOfValues();
 256   CPPUNIT_ASSERT_EQUAL(nbOfValues, theField_2->getNumberOfValues());
 257
 258   if (isFIELD)
 259     {
 260       // Value type and Interlacing type
 261       CPPUNIT_ASSERT_EQUAL(theField_1->getValueType(), theField_2->getValueType());
 262       CPPUNIT_ASSERT_EQUAL(theField_1->getInterlacingType(), theField_2->getInterlacingType());
 263
 264       // Gauss Presence
 265       if (isValue)
 266         {
 267           CPPUNIT_ASSERT_EQUAL(theField_1->getGaussPresence(), theField_2->getGaussPresence());
 268         }
 269       else
 270         {
 271           CPPUNIT_ASSERT_THROW(theField_1->getGaussPresence(), MEDEXCEPTION);
 272           CPPUNIT_ASSERT_THROW(theField_2->getGaussPresence(), MEDEXCEPTION);
 273         }
 274     }
 275   else
 276     {
 277       CPPUNIT_ASSERT_THROW(theField_1->getGaussPresence(), MEDEXCEPTION);
 278       CPPUNIT_ASSERT_THROW(theField_2->getGaussPresence(), MEDEXCEPTION);
 279     }
 280 }
 281
 282 static void checkField_(FIELD_ * theField_, const SUPPORT * theSupport,
 283                         MED_EN::med_type_champ theValueType,
 284                         MED_EN::medModeSwitch theInterlace)
 285 {
 286   // name
 287   const string aFieldName = "a_name_of_a_field";
 288   theField_->setName(aFieldName);
 289   CPPUNIT_ASSERT_EQUAL(aFieldName, theField_->getName());
 290
 291   // description
 292   const string aFieldDescr = "a_description_of_a_field";
 293   theField_->setDescription(aFieldDescr);
 294   CPPUNIT_ASSERT_EQUAL(aFieldDescr, theField_->getDescription());
 295
 296   // support
 297   theField_->setSupport(theSupport);
 298   CPPUNIT_ASSERT(theField_->getSupport() == theSupport);
 299
 300   // components information
 301   int aNbComps = 3;
 302
 303   string aCompsNames[3] =
 304     {
 305       "Vx", "Vy", "Vz"
 306     };
 307   string aCompsDescs[3] =
 308     {
 309       "vitesse selon x", "vitesse selon y", "vitesse selon z"
 310     };
 311   string aCompsUnits[3] =
 312     {
 313       "m.s-1", "m.s-1", "m.s-1"
 314     };
 315
 316   theField_->setNumberOfComponents(aNbComps);
 317   CPPUNIT_ASSERT_EQUAL(aNbComps, theField_->getNumberOfComponents());
 318
 319   theField_->setComponentsNames(aCompsNames);
 320
 321   //#ifdef ENABLE_FAULTS
 322   try
 323     {
 324       theField_->setNumberOfComponents(7);
 325       // Segmentation fault here because array of components names is not resized
 326       for (int i = 1; i <= 7; i++)
 327         {
 328           theField_->setComponentName(i, "AnyComponent");
 329         }
 330     }
 331   catch (MEDEXCEPTION& ex)
 332     {
 333       // Ok, it is good to have MEDEXCEPTION here
 334     }
 335   catch (...)
 336     {
 337       CPPUNIT_FAIL("Unknown exception cought");
 338     }
 339   // restore components names
 340   theField_->setNumberOfComponents(aNbComps);
 341   theField_->setComponentsNames(aCompsNames);
 342   //#endif
 343   //#ifdef ENABLE_FORCED_FAILURES
 344   //  CPPUNIT_FAIL("FIELD_::_componentsNames bad management");
 345   //#endif
 346
 347   theField_->setComponentsDescriptions(aCompsDescs);
 348   theField_->setMEDComponentsUnits(aCompsUnits);
 349
 350   const string * aCompsNamesBack = theField_->getComponentsNames();
 351   const string * aCompsDescsBack = theField_->getComponentsDescriptions();
 352   const string * aCompsUnitsBack = theField_->getMEDComponentsUnits();
 353   for (int i = 1; i <= aNbComps; i++)
 354     {
 355       CPPUNIT_ASSERT_EQUAL(aCompsNamesBack[i-1], theField_->getComponentName(i));
 356       CPPUNIT_ASSERT_EQUAL(aCompsNamesBack[i-1], aCompsNames[i-1]);
 357
 358       CPPUNIT_ASSERT_EQUAL(aCompsDescsBack[i-1], theField_->getComponentDescription(i));
 359       CPPUNIT_ASSERT_EQUAL(aCompsDescsBack[i-1], aCompsDescs[i-1]);
 360
 361       CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack[i-1], theField_->getMEDComponentUnit(i));
 362       CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack[i-1], aCompsUnits[i-1]);
 363     }
 364
 365   const string aCompName2 ("Name of second component");
 366   const string aCompDesc2 ("Description of second component");
 367   const string aCompUnit2 ("Unit of second MED component");
 368
 369   theField_->setComponentName(2, aCompName2);
 370   theField_->setComponentDescription(2, aCompDesc2);
 371   theField_->setMEDComponentUnit(2, aCompUnit2);
 372
 373   const string * aCompsNamesBack2 = theField_->getComponentsNames();
 374   const string * aCompsDescsBack2 = theField_->getComponentsDescriptions();
 375   const string * aCompsUnitsBack2 = theField_->getMEDComponentsUnits();
 376
 377   CPPUNIT_ASSERT_EQUAL(aCompsNamesBack2[1], theField_->getComponentName(2));
 378   CPPUNIT_ASSERT_EQUAL(aCompsNamesBack2[1], aCompName2);
 379
 380   CPPUNIT_ASSERT_EQUAL(aCompsDescsBack2[1], theField_->getComponentDescription(2));
 381   CPPUNIT_ASSERT_EQUAL(aCompsDescsBack2[1], aCompDesc2);
 382
 383   CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack2[1], theField_->getMEDComponentUnit(2));
 384   CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack2[1], aCompUnit2);
 385
 386   //#ifdef ENABLE_FAULTS
 387   // (BUG) No index checking
 388   CPPUNIT_ASSERT_THROW(theField_->setComponentName(0, "str"), MEDEXCEPTION);
 389   CPPUNIT_ASSERT_THROW(theField_->setComponentName(aNbComps + 1, "str"), MEDEXCEPTION);
 390   CPPUNIT_ASSERT_THROW(theField_->setComponentDescription(0, "str"), MEDEXCEPTION);
 391   CPPUNIT_ASSERT_THROW(theField_->setComponentDescription(aNbComps + 1, "str"), MEDEXCEPTION);
 392   CPPUNIT_ASSERT_THROW(theField_->setMEDComponentUnit(0, "str"), MEDEXCEPTION);
 393   CPPUNIT_ASSERT_THROW(theField_->setMEDComponentUnit(aNbComps + 1, "str"), MEDEXCEPTION);
 394   //#endif
 395   //#ifdef ENABLE_FORCED_FAILURES
 396   //  CPPUNIT_FAIL("FIELD::setComponentXXX() does not check component index");
 397   //#endif
 398
 399   // iteration information
 400   int anIterNumber = 10; // set value to MED_NOPDT if undefined (default)
 401   theField_->setIterationNumber(anIterNumber);
 402   CPPUNIT_ASSERT_EQUAL(anIterNumber, theField_->getIterationNumber());
 403
 404   int anOrderNumber = 1; // set value to MED_NONOR if undefined (default)
 405   theField_->setOrderNumber(anOrderNumber);
 406   CPPUNIT_ASSERT_EQUAL(anOrderNumber, theField_->getOrderNumber());
 407
 408   double aTime = 3.435678; // in second
 409   theField_->setTime(aTime);
 410   CPPUNIT_ASSERT_DOUBLES_EQUAL(aTime, theField_->getTime(), 0.0000001);
 411
 412   // Value
 413   int nbOfValues = 10;
 414   // dangerous method, because it does not reallocate values array
 415   theField_->setNumberOfValues(nbOfValues);
 416   CPPUNIT_ASSERT_EQUAL(nbOfValues, theField_->getNumberOfValues());
 417
 418   // Value type and Interlacing type
 419   CPPUNIT_ASSERT_EQUAL(theValueType, theField_->getValueType());
 420   CPPUNIT_ASSERT_EQUAL(theInterlace, theField_->getInterlacingType());
 421 }
 422
 423 template<class T, class INTERLACING_TAG>
 424 void compareField(const FIELD<T, INTERLACING_TAG> * theField_1,
 425                   const FIELD<T, INTERLACING_TAG> * theField_2, bool isValue)
 426 {
 427   // compare FIELD_ part
 428   compareField_(theField_1, theField_2, /*isFIELD = */true, isValue);
 429
 430   // compare FIELD part
 431   // TO DO
 432 }
 433
 434 template<class T, class INTERLACING_TAG>
 435 void checkField (FIELD<T, INTERLACING_TAG> * theField, const SUPPORT * theSupport)
 436 {
 437   // check FIELD_ part
 438   MED_EN::med_type_champ aValueType = SET_VALUE_TYPE<T>::_valueType;
 439   MED_EN::medModeSwitch  anInterlace = SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
 440   checkField_(theField, theSupport, aValueType, anInterlace);
 441
 442   // check FIELD part
 443
 444   // filling by support charackteristics (NOT IMPLEMENTED METHODS!!!):
 445   // value type must be MED_REEL64 (i.e. a FIELD<double>) for these methods,
 446   // nb. of components must be equal 1 (for Volume, Area, Length) or
 447   // space dimension (for Normal, Barycenter, )
 448   {
 449     const GMESH* aMesh = theSupport->getMesh();
 450     int spaceDim = 3;
 451     if (aMesh) spaceDim = aMesh->getSpaceDimension();
 452     theField->deallocValue();
 453     theField->allocValue(/*NumberOfComponents = */spaceDim + 1);
 454
 455     //  0020142: [CEA 315] Unused function in MEDMEM::FIELD
 456     // getVolume() etc. does nothing
 457     //     CPPUNIT_ASSERT_THROW(theField->getVolume(), MEDEXCEPTION);
 458     //     CPPUNIT_ASSERT_THROW(theField->getArea(), MEDEXCEPTION);
 459     //     CPPUNIT_ASSERT_THROW(theField->getLength(), MEDEXCEPTION);
 460     //     if (aMesh) {
 461     //       CPPUNIT_ASSERT_THROW(theField->getNormal(), MEDEXCEPTION);
 462     //       CPPUNIT_ASSERT_THROW(theField->getBarycenter(), MEDEXCEPTION);
 463     //     }
 464
 465     theField->deallocValue();
 466     theField->allocValue(/*NumberOfComponents = */1);
 467     //  0020142: [CEA 315] Unused function in MEDMEM::FIELD
 468     // getVolume() etc. does nothing
 469     //     if (aValueType == MED_EN::MED_REEL64) {
 470     //       CPPUNIT_ASSERT_NO_THROW(theField->getVolume());
 471     //       CPPUNIT_ASSERT_NO_THROW(theField->getArea());
 472     //       CPPUNIT_ASSERT_NO_THROW(theField->getLength());
 473     //     }
 474     //     else {
 475     //       CPPUNIT_ASSERT_THROW(theField->getVolume(), MEDEXCEPTION);
 476     //       CPPUNIT_ASSERT_THROW(theField->getArea(), MEDEXCEPTION);
 477     //       CPPUNIT_ASSERT_THROW(theField->getLength(), MEDEXCEPTION);
 478     //     }
 479
 480     if (aMesh)
 481       {
 482         theField->deallocValue();
 483         theField->allocValue(/*NumberOfComponents = */spaceDim);
 484         //  0020142: [CEA 315] Unused function in MEDMEM::FIELD
 485         // getVolume() etc. does nothing
 486         //       if (aValueType == MED_EN::MED_REEL64) {
 487         //         CPPUNIT_ASSERT_NO_THROW(theField->getNormal());
 488         //         CPPUNIT_ASSERT_NO_THROW(theField->getBarycenter());
 489         //       }
 490         //       else {
 491         //         CPPUNIT_ASSERT_THROW(theField->getNormal(), MEDEXCEPTION);
 492         //         CPPUNIT_ASSERT_THROW(theField->getBarycenter(), MEDEXCEPTION);
 493         //       }
 494       }
 495   }
 496
 497   // values
 498   theField->deallocValue();
 499   theField->allocValue(/*NumberOfComponents = */2);
 500   int nbElemSupport = theSupport->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
 501   CPPUNIT_ASSERT_EQUAL(nbElemSupport, theField->getNumberOfValues());
 502
 503   //#ifdef ENABLE_FAULTS
 504   // (BUG) FIELD::deallocValue() does not nullify _value pointer,
 505   // that is why there can be failures in other methods
 506   // (even if simply call deallocValue() two times)
 507   theField->deallocValue();
 508   theField->getGaussPresence();
 509   //#endif
 510   //#ifdef ENABLE_FORCED_FAILURES
 511   //  CPPUNIT_FAIL("FIELD::deallocValue() does not nullify _value pointer");
 512   //#endif
 513
 514   // copy constructor
 515   FIELD<T, INTERLACING_TAG> aField_copy1 (*theField);
 516   //compareField(theField, &aField_copy1, /*isValue = */false);
 517   compareField(theField, &aField_copy1, /*isValue = */true);
 518
 519   // operator=
 520   //#ifdef ENABLE_FAULTS
 521   // (BUG) This fails (Segmentation fault) if not set:
 522   // _componentsNames or _componentsDescriptions, or _componentsUnits, or _MEDComponentsUnits
 523   FIELD<T, INTERLACING_TAG> aField_copy2;
 524   aField_copy2 = *theField;
 525   //compareField(theField, &aField_copy2, /*isValue = */false);
 526   compareField(theField, &aField_copy2, /*isValue = */true);
 527   //#endif
 528   //#ifdef ENABLE_FORCED_FAILURES
 529   //  CPPUNIT_FAIL("FIELD_::operator=() fails if _componentsUnits is not set");
 530   //#endif
 531 }
 532
 533 template<class T>
 534 FIELD<T> * createFieldOnGroup(MESH* theMesh, const GROUP* theGroup,
 535                               const string theName, const string theDescr)
 536 {
 537   FIELD<T> * aFieldOnGroup = new FIELD<T> (theGroup, /*NumberOfComponents = */2);
 538
 539   aFieldOnGroup->setName(theName);
 540   aFieldOnGroup->setDescription(theDescr);
 541
 542   string aCompsNames[2] =
 543     {
 544       "Pos", "Neg"
 545     };
 546   string aCompsDescs[2] =
 547     {
 548       "+", "-"
 549     };
 550   string aCompsUnits[2] =
 551     {
 552       "unit1", "unit2"
 553     };
 554
 555   aFieldOnGroup->setComponentsNames(aCompsNames);
 556   aFieldOnGroup->setComponentsDescriptions(aCompsDescs);
 557   aFieldOnGroup->setMEDComponentsUnits(aCompsUnits);
 558
 559   return aFieldOnGroup;
 560 }
 561
 562 double plus13 (double val);
 563 double plus13 (double val)
 564 {
 565   return val + 13;
 566 }
 567
 568 // function to calculate field values from coordinates of an element
 569 // typedef void (*myFuncType)(const double * temp, T* output);
 570 // size of temp array = space dim = 3
 571 // size of output array = nb. comps = 2
 572 static void proj2d (const double * temp, double* output)
 573 {
 574   // dimetric projection with coefficients:
 575   // 1.0 along Oy and Oz, 0.5 along Ox
 576   //
 577   //    ^ z (y_)
 578   //    |
 579   //    |
 580   //    .----> y (x_)
 581   //   /
 582   //  L x
 583   //
 584   // x_ = y - x * sqrt(2.) / 4.
 585   // y_ = z - x * sqrt(2.) / 4.
 586
 587   double dx = temp[0] * std::sqrt(2.) / 4.;
 588   output[0] = temp[1] - dx;
 589   output[1] = temp[2] - dx;
 590 }
 591
 592 static void testDrivers()
 593 {
 594   string filename_rd                  = getResourceFile("pointe.med");
 595   string filename_wr                  = makeTmpFile("myMedFieldfile.med", filename_rd);
 596   string filename_support_wr          = makeTmpFile("myMedSupportFiledfile.med");
 597   string filename22_rd                = getResourceFile("pointe.med");
 598   string filenamevtk_wr               = makeTmpFile("myMedFieldfile22.vtk");
 599
 600   string fieldname_celldouble_rd      = "fieldcelldoublescalar";
 601   string fieldname_celldouble_wr      = fieldname_celldouble_rd + "_cpy";
 602   string fieldname_nodeint_rd         = "fieldnodeint";
 603   string fieldname_nodeint_wr         = fieldname_nodeint_rd + "_cpy";
 604   string fieldname_nodeint_wr1        = fieldname_nodeint_rd + "_cpy1";
 605   string meshname                     = "maa1";
 606
 607   // To remove tmp files from disk
 608   MEDMEMTest_TmpFilesRemover aRemover;
 609   aRemover.Register(filename_wr);
 610   aRemover.Register(filenamevtk_wr);
 611   aRemover.Register(filename_support_wr);
 612
 613   FIELD<int> *aInvalidField=new FIELD<int>();
 614   //must throw becase only VTK_DRIVER or MED_DRIVER may be specified as driverType for FIELD
 615   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(NO_DRIVER, filename_rd, fieldname_nodeint_rd)),
 616                        MEDEXCEPTION);
 617   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(GIBI_DRIVER, filename_rd, fieldname_nodeint_rd)),
 618                        MEDEXCEPTION);
 619   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(PORFLOW_DRIVER, filename_rd, fieldname_nodeint_rd)),
 620                        MEDEXCEPTION);
 621   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(ASCII_DRIVER, filename_rd, fieldname_nodeint_rd)),
 622                        MEDEXCEPTION);
 623
 624   //////////////////
 625   //TestRead Part//
 626   //////////////////
 627   FIELD<double> *aField_1 = NULL;
 628   CPPUNIT_ASSERT_NO_THROW(aField_1 = new FIELD<double>(MED_DRIVER, filename_rd, fieldname_celldouble_rd));
 629
 630   //Test read(int index) method
 631   int IdDriver_rd = aField_1->addDriver(MED_DRIVER,filename_rd,fieldname_celldouble_rd);
 632   //#ifdef ENABLE_FORCED_FAILURES
 633   // (BUG) Cannot open file, but file exist
 634   CPPUNIT_ASSERT_NO_THROW(aField_1->read(IdDriver_rd));
 635   //#endif
 636
 637   //Test read(GENDRIVER & genDriver) method
 638   //Creation a Driver
 639   MED_FIELD_RDONLY_DRIVER<int> *aMedRdFieldDriver_1 = new MED_FIELD_RDONLY_DRIVER<int>();
 640   //Creation a Field
 641   FIELD<int> *aField_2 = new FIELD<int>();
 642   aField_2->setName(fieldname_nodeint_rd);
 643   aField_2->addDriver(*aMedRdFieldDriver_1);
 644   aField_2->read(*aMedRdFieldDriver_1);
 645
 646   ///////////////////
 647   //Test Write Part//
 648   ///////////////////
 649   int IdDriver;
 650   MESH *aMesh = new MESH(MED_DRIVER,filename_rd,meshname);
 651   const SUPPORT *aSupport = aMesh->getSupportOnAll( MED_CELL );
 652   FIELD<int> *aFieldSupport;
 653   //#ifdef ENABLE_FORCED_FAILURES
 654   CPPUNIT_ASSERT_NO_THROW(aFieldSupport =
 655                           new FIELD<int>(aSupport, MED_DRIVER,filename_support_wr,fieldname_nodeint_rd));
 656   //(BUG) Can not open file
 657   MED_FIELD_WRONLY_DRIVER<int> * aFieldWrDriver =
 658     new MED_FIELD_WRONLY_DRIVER<int>(filename_support_wr,aFieldSupport);
 659   aFieldWrDriver->setFieldName(aFieldSupport->getName() + "_copy");
 660   CPPUNIT_ASSERT_NO_THROW(IdDriver= aFieldSupport->addDriver(*aFieldWrDriver));
 661   CPPUNIT_ASSERT_NO_THROW(aFieldSupport->write(IdDriver));
 662   aFieldSupport->removeReference();
 663   delete aFieldWrDriver;
 664   //#endif
 665
 666   //Create fileds
 667   FIELD<double> * aField_3 = new FIELD<double>();
 668   MED_FIELD_RDONLY_DRIVER<double> *aMedRdFieldDriver_2 =
 669     new MED_FIELD_RDONLY_DRIVER<double>(filename_rd, aField_3);
 670   aMedRdFieldDriver_2->open();
 671   aMedRdFieldDriver_2->setFieldName(fieldname_celldouble_rd);
 672   aMedRdFieldDriver_2->read();
 673   aMedRdFieldDriver_2->close();
 674
 675   //Test write(int index) method
 676   //Add drivers to FIELDs
 677   int IdDriver1 = -1;
 678   try
 679     {
 680       IdDriver1 = aField_3->addDriver(MED_DRIVER,filename_wr,fieldname_celldouble_wr);
 681     }
 682   catch(MEDEXCEPTION &e)
 683     {
 684       e.what();
 685     }
 686   catch( ... )
 687     {
 688       CPPUNIT_FAIL("Unknown exception");
 689     }
 690   //Trying call write(int index) method with incorrect index
 691   //#ifdef ENABLE_FAULTS
 692   CPPUNIT_ASSERT_THROW(aField_3->write(IdDriver1+1),MEDEXCEPTION);
 693   // => Segmentation fault
 694   //#endif
 695
 696   //Write field to file
 697   //#ifdef ENABLE_FAULTS
 698   try
 699     {
 700       aField_3->write(IdDriver1);
 701       // => Segmentation fault
 702     }
 703   catch(MEDEXCEPTION &e)
 704     {
 705       e.what();
 706     }
 707   catch( ... )
 708     {
 709       CPPUNIT_FAIL("Unknown exception");
 710     }
 711   //#endif
 712
 713   CPPUNIT_ASSERT_NO_THROW(aField_3->rmDriver(IdDriver1));
 714
 715   //Test write(const GENDRIVER &);
 716   //Create a driver
 717   MED_FIELD_WRONLY_DRIVER<int> *aMedWrFieldDriver = new MED_FIELD_WRONLY_DRIVER<int>();
 718   aMedWrFieldDriver->setFileName(filename_wr);
 719   aField_2->setName(fieldname_nodeint_wr1);
 720   //Add driver to a field
 721   aField_2->addDriver(*aMedWrFieldDriver);
 722
 723   try
 724     {
 725       aField_2->write(*aMedWrFieldDriver);
 726     }
 727   catch(MEDEXCEPTION &e)
 728     {
 729       e.what();
 730     }
 731   catch( ... )
 732     {
 733       CPPUNIT_FAIL("Unknown exception");
 734     }
 735
 736   //Test writeAppend(int index) method
 737   //Create a vtk file
 738   MESH * aMesh_1 = new MESH;
 739   MED_MESH_RDONLY_DRIVER *aMedMeshRdDriver22 = new MED_MESH_RDONLY_DRIVER(filename22_rd, aMesh_1);
 740   aMedMeshRdDriver22->open();
 741   aMedMeshRdDriver22->setMeshName(meshname);
 742   aMedMeshRdDriver22->read();
 743   aMedMeshRdDriver22->close();
 744   VTK_MESH_DRIVER *aVtkDriver = new VTK_MESH_DRIVER(filenamevtk_wr, aMesh_1);
 745   aVtkDriver->open();
 746   aVtkDriver->write();
 747   aVtkDriver->close();
 748
 749   //Create a field
 750   FIELD<int> * aField_4 = new FIELD<int>();
 751   MED_FIELD_RDONLY_DRIVER<int> *aMedRdFieldDriver22 =
 752     new MED_FIELD_RDONLY_DRIVER<int>(filename22_rd, aField_2);
 753   aMedRdFieldDriver22->open();
 754   aMedRdFieldDriver22->setFieldName(fieldname_nodeint_rd);
 755   aMedRdFieldDriver22->read();
 756   aMedRdFieldDriver22->close();
 757
 758   //Add Driver to a field
 759   int IdDriver2;
 760   try
 761     {
 762       IdDriver2 = aField_4->addDriver(VTK_DRIVER, filenamevtk_wr ,fieldname_nodeint_wr);
 763     }
 764   catch(MEDEXCEPTION &e)
 765     {
 766       e.what();
 767     }
 768   catch( ... )
 769     {
 770       CPPUNIT_FAIL("Unknown exception");
 771     }
 772   //#ifdef ENABLE_FAULTS
 773   //Trying call writeAppend() method with incorrect index
 774   CPPUNIT_ASSERT_THROW(aField_4->writeAppend(IdDriver2+1,fieldname_nodeint_wr),MEDEXCEPTION);
 775   // => Segmentation fault
 776   //#endif
 777
 778   //#ifdef ENABLE_FAULTS
 779   // (BUG) => Segmentation fault
 780   CPPUNIT_ASSERT_NO_THROW(aField_4->writeAppend(IdDriver2, fieldname_nodeint_wr));
 781   //#endif
 782
 783   //Test writeAppend(const GENDRIVER &) method
 784   aField_4->setName(fieldname_nodeint_wr1);
 785
 786   //Add driver to a field
 787   //#ifdef ENABLE_FAULTS
 788   //Create a driver
 789   VTK_FIELD_DRIVER<int> *aVtkFieldDriver = new VTK_FIELD_DRIVER<int>(filenamevtk_wr, aField_4);
 790   CPPUNIT_ASSERT_NO_THROW(aField_4->addDriver(*aVtkFieldDriver));
 791   //(BUG) => Segmentation fault after addDriver(const GENDRIVER &)
 792   CPPUNIT_ASSERT_NO_THROW(aField_4->writeAppend(*aVtkFieldDriver));
 793   //#endif
 794
 795
 796   //Delete objects
 797   aField_1->removeReference();
 798   delete aMedRdFieldDriver_1;
 799   aField_2->removeReference();
 800   aField_3->removeReference();
 801   delete aMedRdFieldDriver_2;
 802   aField_4->removeReference();
 803   delete aMedMeshRdDriver22;
 804   delete aMedWrFieldDriver;
 805   delete aVtkDriver;
 806   aMesh->removeReference();
 807   aMesh_1->removeReference();
 808   delete aMedRdFieldDriver22;
 809 }
 810
 811 static void MEDMEMTest_testField()
 812 {
 813   SUPPORT *anEmptySupport=new SUPPORT;
 814   ////////////////////
 815   // TEST 1: FIELD_ //
 816   ////////////////////
 817   FIELD_ *aField_=new FIELD_ ;
 818
 819   // check set/get methods
 820   MED_EN::med_type_champ aValueType = MED_EN::MED_UNDEFINED_TYPE;
 821   MED_EN::medModeSwitch  anInterlace = MED_EN::MED_UNDEFINED_INTERLACE;
 822   checkField_(aField_, anEmptySupport, aValueType, anInterlace);
 823
 824   // copy constructor
 825   // This fails (Segmentation fault) if not set:
 826   // _componentsNames or _componentsDescriptions, or _MEDComponentsUnits
 827   FIELD_ *aField_copy1=new FIELD_(*aField_);
 828   compareField_(aField_, aField_copy1, /*isFIELD = */false, /*isValue = */false);
 829
 830   // operator=
 831   //#ifdef ENABLE_FAULTS
 832   // (BUG) This fails (Segmentation fault) if not set:
 833   // _componentsNames or _componentsDescriptions, or _componentsUnits, or _MEDComponentsUnits
 834   // (BUG) Code duplication with copyGlobalInfo(), called from copy constructor
 835   FIELD_ *aField_copy2=new FIELD_;
 836   *aField_copy2 = *aField_;
 837   compareField_(aField_, aField_copy2, /*isFIELD = */false, /*isValue = */false);
 838   //#endif
 839   //#ifdef ENABLE_FORCED_FAILURES
 840   //  CPPUNIT_FAIL("FIELD_::operator=() fails if _componentsUnits is not set");
 841   //#endif
 842
 843   // following test is commented since method
 844   // setTotalNumberOfElements() is removed.
 845   /*
 846   // construction on a given support
 847   {
 848   anEmptySupport.setTotalNumberOfElements(11);
 849   // CASE1:
 850   FIELD_ aField_case1 (&anEmptySupport, 10);
 851   // CASE2:
 852   FIELD_ aField_case2;
 853   aField_case2.setSupport(&anEmptySupport);
 854   aField_case2.setNumberOfComponents(10);
 855
 856   #ifdef ENABLE_FORCED_FAILURES
 857   CPPUNIT_ASSERT_EQUAL_MESSAGE("No correspondance between CASE1 and CASE2",
 858   aField_case1.getNumberOfValues(),
 859   aField_case2.getNumberOfValues());
 860   #endif
 861   }
 862   */
 863
 864   ////////////////////////
 865   // TEST 2: FIELD<int> //
 866   ////////////////////////
 867   FIELD<int> *aFieldInt=new FIELD<int>();
 868   checkField(aFieldInt, anEmptySupport);
 869
 870   ////////////////////////////////////////
 871   // TEST 3: FIELD<double, NoInterlace> //
 872   ////////////////////////////////////////
 873   MESH * aMesh = MEDMEMTest_createTestMesh();
 874   const GROUP* aGroup = aMesh->getGroup(MED_EN::MED_FACE, 1);
 875
 876   FIELD<double, NoInterlace> aFieldDouble;
 877   checkField(&aFieldDouble, aGroup);
 878
 879   //////////////////////////////////////////
 880   // TEST 4: FIELD<double, FullInterlace> //
 881   //////////////////////////////////////////
 882   FIELD<double> * aFieldOnGroup1 = createFieldOnGroup<double>(aMesh, aGroup, "Linear", "N");
 883   FIELD<double> * aFieldOnGroup2 = createFieldOnGroup<double>(aMesh, aGroup, "Quadratic", "N**2");
 884
 885   int nbVals = aFieldOnGroup1->getNumberOfValues();
 886   CPPUNIT_ASSERT(nbVals);
 887
 888   // numbers of elements in group,
 889   // they are needed in method FIELD::setValueIJ()
 890   const int *anElems = aGroup->getnumber()->getValue();
 891   double eucl1 = 0., eucl2 = 0.;
 892
 893   for (int i = 1; i <= nbVals; i++)
 894     {
 895       aFieldOnGroup1->setValueIJ(anElems[i-1], 1, (double)i);
 896       aFieldOnGroup1->setValueIJ(anElems[i-1], 2, (double)(-i));
 897
 898       aFieldOnGroup2->setValueIJ(anElems[i-1], 1, (double)i*i);
 899       aFieldOnGroup2->setValueIJ(anElems[i-1], 2, (double)(-i*i));
 900
 901       eucl1 += 2. * i * i;
 902       eucl2 += 2. * i * i * i * i;
 903     }
 904
 905   // out of bound (inexisting 33-th component of last element)
 906   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->setValueIJ(anElems[nbVals-1], 33, 10.), MEDEXCEPTION);
 907
 908   // normMax
 909   CPPUNIT_ASSERT_DOUBLES_EQUAL(nbVals, aFieldOnGroup1->normMax(), 0.000001);
 910   CPPUNIT_ASSERT_DOUBLES_EQUAL(nbVals*nbVals, aFieldOnGroup2->normMax(), 0.000001);
 911
 912   // norm2
 913   CPPUNIT_ASSERT_DOUBLES_EQUAL(std::sqrt(eucl1), aFieldOnGroup1->norm2(), 0.000001); // 10.4881
 914   CPPUNIT_ASSERT_DOUBLES_EQUAL(std::sqrt(eucl2), aFieldOnGroup2->norm2(), 0.000001); // 44.2493
 915
 916   // check getXXX methods
 917   CPPUNIT_ASSERT(!aFieldOnGroup1->getGaussPresence());
 918   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getArrayGauss(), MEDEXCEPTION);
 919   MEDMEM_ArrayInterface<double, FullInterlace,NoGauss>::Array * anArrayNoGauss =
 920     aFieldOnGroup1->getArrayNoGauss();
 921
 922   MEDMEM_Array_ * aMEDMEM_Array_ = aFieldOnGroup1->getArray();
 923   MEDMEM_ArrayInterface<double, FullInterlace,NoGauss>::Array * aMEDMEM_Array_conv =
 924     static_cast<MEDMEM_ArrayInterface<double, FullInterlace,NoGauss>::Array *>(aMEDMEM_Array_);
 925
 926   const double * aValues = aFieldOnGroup1->getValue();
 927
 928   // out of range
 929   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getColumn(3), MEDEXCEPTION);
 930   // cannot get column in FullInterlace
 931   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getColumn(1), MEDEXCEPTION);
 932
 933   for (int i = 1; i <= nbVals; i++)
 934     {
 935       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , aFieldOnGroup1->getValueIJK(anElems[i-1], 1, 1), 0.000001);
 936       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), aFieldOnGroup1->getValueIJK(anElems[i-1], 2, 1), 0.000001);
 937
 938       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , aValues[(i-1)*2 + 0], 0.000001);
 939       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), aValues[(i-1)*2 + 1], 0.000001);
 940
 941       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , anArrayNoGauss->getIJ(i, 1), 0.000001);
 942       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), anArrayNoGauss->getIJ(i, 2), 0.000001);
 943
 944       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , aMEDMEM_Array_conv->getIJ(i, 1), 0.000001);
 945       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), aMEDMEM_Array_conv->getIJ(i, 2), 0.000001);
 946
 947       const double* row_i = aFieldOnGroup1->getRow(anElems[i-1]);
 948       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , row_i[0], 0.000001);
 949       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), row_i[1], 0.000001);
 950
 951       double vals_i [2];
 952       aFieldOnGroup1->getValueOnElement(anElems[i-1], vals_i);
 953       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , vals_i[0], 0.000001);
 954       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), vals_i[1], 0.000001);
 955     }
 956
 957   // modify all values of aFieldOnGroup2 by formula a*x + b (a = 2, b = 3)
 958   aFieldOnGroup2->applyLin(2., 3.);
 959   for (int i = 1; i <= nbVals; i++)
 960     {
 961       CPPUNIT_ASSERT_DOUBLES_EQUAL(3. + 2.*i*i, aFieldOnGroup2->getValueIJ(anElems[i-1], 1), 0.000001);
 962       CPPUNIT_ASSERT_DOUBLES_EQUAL(3. - 2.*i*i, aFieldOnGroup2->getValueIJ(anElems[i-1], 2), 0.000001);
 963     }
 964
 965   // apply function plus13() to aFieldOnGroup1
 966   aFieldOnGroup1->applyFunc<plus13>();
 967   for (int i = 1; i <= nbVals; i++)
 968     {
 969       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
 970       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
 971     }
 972
 973   // scalarProduct
 974   FIELD<double, FullInterlace> * aScalarProduct =
 975     FIELD<double, FullInterlace>::scalarProduct(*aFieldOnGroup1, *aFieldOnGroup2, /*deepCheck = */true);
 976   CPPUNIT_ASSERT_EQUAL(nbVals, aScalarProduct->getNumberOfValues());
 977   CPPUNIT_ASSERT_EQUAL(1, aScalarProduct->getNumberOfComponents());
 978   for (int i = 1; i <= nbVals; i++)
 979     {
 980       CPPUNIT_ASSERT_DOUBLES_EQUAL(78. + 4.*i*i*i, //(3. + 2.*i*i)*(13 + i) + (3. - 2.*i*i)*(13 - i)
 981                                    aScalarProduct->getValueIJ(anElems[i-1], 1), 0.000001);
 982     }
 983
 984   // fillFromAnalytic
 985   aFieldOnGroup2->fillFromAnalytic(proj2d);
 986
 987   double bary [3];
 988   double outp [2];
 989   const SUPPORT * aSupp = aFieldOnGroup2->getSupport();
 990   FIELD<double,FullInterlace> * barycenter = aMesh->getBarycenter(aSupp);
 991   for (int i = 1; i <= nbVals; i++)
 992     {
 993       bary[0] = barycenter->getValueIJ(anElems[i-1], 1);
 994       bary[1] = barycenter->getValueIJ(anElems[i-1], 2);
 995       bary[2] = barycenter->getValueIJ(anElems[i-1], 3);
 996
 997       proj2d(bary, outp);
 998
 999       //cout << "barycenter (" << bary[0] << ", " << bary[1] << ", " << bary[2] << ")" << endl;
1000       //cout << "proj2d     (" << outp[0] << ", " << outp[1] << ")" << endl;
1001
1002       //bary (-0.666667,  0.666667, 0.666667) -> outp ( 0.902369, 0.902369)
1003       //bary ( 0.666667, -0.666667, 0.666667) -> outp (-0.902369, 0.430964)
1004       //bary ( 0.      ,  0.      , 2.      ) -> outp ( 0.      , 2.      )
1005       //bary ( 0.      ,  0.      , 3.      ) -> outp ( 0.      , 3.      )
1006       //bary (-1.      ,  0.      , 2.5     ) -> outp ( 0.353553, 2.85355 )
1007
1008       //#ifdef ENABLE_FORCED_FAILURES
1009       // (BUG) in FIELD::fillFromAnalytic() in case of support, different from nodes:
1010       //       barycenterField in FullInterlace, but values extracted like from NoInterlace
1011       CPPUNIT_ASSERT_DOUBLES_EQUAL(outp[0], aFieldOnGroup2->getValueIJ(anElems[i-1], 1), 0.000001);
1012       CPPUNIT_ASSERT_DOUBLES_EQUAL(outp[1], aFieldOnGroup2->getValueIJ(anElems[i-1], 2), 0.000001);
1013       //#endif
1014
1015       // currently it gives values, that are wrong:
1016       //aFieldOnGroup2 row1 ( 0.902369,  0.235702)
1017       //aFieldOnGroup2 row2 (-0.235702,  2.7643  )
1018       //aFieldOnGroup2 row3 (-0.235702, -1.2357  )
1019       //aFieldOnGroup2 row4 ( 1.7643  , -0.235702)
1020       //aFieldOnGroup2 row5 ( 0.235702,  2.7357  )
1021     }
1022
1023   // info about support (Group1)
1024   CPPUNIT_ASSERT(!aFieldOnGroup1->isOnAllElements()); // because we build Group1 so
1025   int nbTypes = aFieldOnGroup1->getNumberOfGeometricTypes();
1026   //CPPUNIT_ASSERT(nbTypes);
1027   CPPUNIT_ASSERT_EQUAL(2, nbTypes);
1028   const int * nbElemsInEachType = aFieldOnGroup1->getNumberOfElements();
1029   const MED_EN::medGeometryElement * aGeomTypes = aFieldOnGroup1->getGeometricTypes();
1030
1031   CPPUNIT_ASSERT_EQUAL(MED_EN::MED_TRIA3, aGeomTypes[0]);
1032   CPPUNIT_ASSERT_EQUAL(MED_EN::MED_QUAD4, aGeomTypes[1]);
1033
1034   // GAUSS
1035
1036   // now we have no gauss localization in aFieldOnGroup1
1037   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA3));
1038   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_QUAD4));
1039   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA6), MEDEXCEPTION);
1040   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(), MEDEXCEPTION);
1041
1042   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalization(MED_EN::MED_TRIA3), MEDEXCEPTION);
1043   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalizationPtr(MED_EN::MED_TRIA3), MEDEXCEPTION);
1044
1045   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNbGaussI(anElems[0]));
1046
1047   // set a gauss localization into aFieldOnGroup1
1048   double cooRef[6] =
1049     {
1050       1.,1., 2.,4., 3.,9.
1051     }; // xy xy xy
1052   double cooGauss[10] =
1053     {
1054       7.,7., 6.,6., 5.,5., 4.,3., 2.,1.
1055     }; // x1,y1  x2,y2  x3,y3  x4,y4  x5,y5
1056   double wg[5] =
1057     {
1058       1., 2., 3., 4., 5.
1059     };
1060   GAUSS_LOCALIZATION<> gl1 ("GL1", MED_EN::MED_TRIA3, /*nGauss*/5, cooRef, cooGauss, wg);
1061
1062   aFieldOnGroup1->setGaussLocalization(MED_EN::MED_TRIA3, gl1);
1063
1064   // now we have a gauss localization for MED_TRIA3 type
1065   CPPUNIT_ASSERT_EQUAL(5, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA3));
1066   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_QUAD4));
1067   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA6), MEDEXCEPTION);
1068   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(), MEDEXCEPTION);
1069
1070   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalization(MED_EN::MED_QUAD4), MEDEXCEPTION);
1071   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalizationPtr(MED_EN::MED_QUAD4), MEDEXCEPTION);
1072
1073   GAUSS_LOCALIZATION<> gl1Back = aFieldOnGroup1->getGaussLocalization(MED_EN::MED_TRIA3);
1074   const GAUSS_LOCALIZATION<> * gl1BackPtr = aFieldOnGroup1->getGaussLocalizationPtr(MED_EN::MED_TRIA3);
1075
1076   CPPUNIT_ASSERT(gl1 == gl1Back);
1077   CPPUNIT_ASSERT(gl1 == *gl1BackPtr);
1078
1079   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNbGaussI(anElems[0]));
1080
1081   // sub-support of Group1 on one (first) geometric type
1082   SUPPORT * aSubSupport1 = new SUPPORT;
1083   aSubSupport1->setMesh( aMesh );
1084   aSubSupport1->setName( "Sub-Support 1 of Group1" );
1085   aSubSupport1->setEntity( MED_EN::MED_FACE );
1086
1087   int nbTypes1 = 1;
1088   int nbElemsInEachType1[1];
1089   nbElemsInEachType1[0] = nbElemsInEachType[0];
1090   int nbElems1 = nbElemsInEachType1[0];
1091   MED_EN::medGeometryElement aGeomTypes1[1];
1092   aGeomTypes1[0] = aGeomTypes[0];
1093   int * anElems1 = new int[nbElems1];
1094   for (int i = 0; i < nbElems1; i++)
1095     {
1096       anElems1[i] = anElems[i];
1097     }
1098
1099   aSubSupport1->setpartial("Support for sub-field 1 on one type of elements",
1100                            nbTypes1, nbElems1, aGeomTypes1, nbElemsInEachType1, anElems1);
1101
1102   //cout << "aSubSupport1:" << endl;
1103   //cout << *aSubSupport1 << endl;
1104
1105   // extract sub-field on aSubSupport1
1106   FIELD<double, FullInterlace> * aSubField1 = aFieldOnGroup1->extract(aSubSupport1);
1107   CPPUNIT_ASSERT_EQUAL(nbElems1 * /*NumberOfComponents = */2, aSubField1->getValueLength());
1108
1109   // aSubField1:
1110   // elt\comp |  1 |  2
1111   //--------------------
1112   //  1       | 14 | 12
1113   //  2       | 15 | 11
1114
1115   // check normL2() and normL1()
1116   FIELD<double>* anAreaField = aMesh->getArea(aSubSupport1);
1117   double area1 = anAreaField->getValueIJ(anElems1[0], 1);
1118   double area2 = anAreaField->getValueIJ(anElems1[1], 1);
1119   CPPUNIT_ASSERT_DOUBLES_EQUAL(2.44949, area1, 0.00001);
1120   CPPUNIT_ASSERT_DOUBLES_EQUAL(2.44949, area2, 0.00001);
1121
1122   CPPUNIT_ASSERT_DOUBLES_EQUAL(210.5, aSubField1->normL2(1), 0.00001); // (14*14 + 15*15)/2
1123   //#ifdef ENABLE_FORCED_FAILURES
1124   // (BUG) FIELD::normL2(int component, const FIELD * p_field_volume):
1125   //       component is not taken into account
1126   CPPUNIT_ASSERT_DOUBLES_EQUAL(132.5, aSubField1->normL2(2), 0.00001); // (12*12 + 11*11)/2
1127   //#endif
1128   CPPUNIT_ASSERT_DOUBLES_EQUAL(343.0, aSubField1->normL2() , 0.00001); // 210.5 + 132.5
1129
1130   CPPUNIT_ASSERT_DOUBLES_EQUAL(14.5, aSubField1->normL1(1), 0.00001); // (14 + 15)/2
1131   CPPUNIT_ASSERT_DOUBLES_EQUAL(11.5, aSubField1->normL1(2), 0.00001); // (12 + 11)/2
1132   CPPUNIT_ASSERT_DOUBLES_EQUAL(26.0, aSubField1->normL1() , 0.00001); // 14.5 + 11.5
1133
1134   double aNewArea [2] =
1135     {
1136       1., 0.
1137     }; // only first element will be taken into account
1138   anAreaField->setValue(aNewArea);
1139
1140   CPPUNIT_ASSERT_DOUBLES_EQUAL(196.0, aSubField1->normL2(1, anAreaField), 0.00001); // 14*14/1
1141   //#ifdef ENABLE_FORCED_FAILURES
1142   // (BUG) FIELD::normL2(int component, const FIELD * p_field_volume):
1143   //       component is not taken into account
1144   CPPUNIT_ASSERT_DOUBLES_EQUAL(144.0, aSubField1->normL2(2, anAreaField), 0.00001); // 12*12/1
1145   //#endif
1146   CPPUNIT_ASSERT_DOUBLES_EQUAL(340.0, aSubField1->normL2(anAreaField) , 0.00001); // 196 + 144
1147
1148   CPPUNIT_ASSERT_DOUBLES_EQUAL(14.0, aSubField1->normL1(1, anAreaField), 0.00001); // 14/1
1149   CPPUNIT_ASSERT_DOUBLES_EQUAL(12.0, aSubField1->normL1(2, anAreaField), 0.00001); // 12/1
1150   CPPUNIT_ASSERT_DOUBLES_EQUAL(26.0, aSubField1->normL1(anAreaField) , 0.00001); // 14 + 12
1151
1152   // applyPow
1153   aSubField1->applyPow(2.);
1154   CPPUNIT_ASSERT_DOUBLES_EQUAL(196., aSubField1->getValueIJ(anElems1[0], 1), 0.000001); // 14*14
1155   CPPUNIT_ASSERT_DOUBLES_EQUAL(144., aSubField1->getValueIJ(anElems1[0], 2), 0.000001); // 12*12
1156   CPPUNIT_ASSERT_DOUBLES_EQUAL(225., aSubField1->getValueIJ(anElems1[1], 1), 0.000001); // 15*15
1157   CPPUNIT_ASSERT_DOUBLES_EQUAL(121., aSubField1->getValueIJ(anElems1[1], 2), 0.000001); // 11*11
1158
1159   // setArray (NoGauss)
1160   MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array * aNewArrayNoGauss =
1161     new MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array(/*dim*/2, /*nbelem*/2);
1162   aNewArrayNoGauss->setIJ(1, 1, 4.);
1163   aNewArrayNoGauss->setIJ(1, 2, 2.);
1164   aNewArrayNoGauss->setIJ(2, 1, 5.);
1165   aNewArrayNoGauss->setIJ(2, 2, 1.);
1166   aSubField1->setArray(aNewArrayNoGauss);
1167   // no need to delete aNewArrayNoGauss, because it will be deleted
1168   // in destructor or in deallocValue() method of aSubField1
1169
1170   CPPUNIT_ASSERT_DOUBLES_EQUAL(4., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
1171   CPPUNIT_ASSERT_DOUBLES_EQUAL(2., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
1172   CPPUNIT_ASSERT_DOUBLES_EQUAL(5., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
1173   CPPUNIT_ASSERT_DOUBLES_EQUAL(1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
1174
1175   // setRow
1176   double row[2] =
1177     {
1178       -1., -3.
1179     };
1180   aSubField1->setRow(anElems1[0], row);
1181   CPPUNIT_ASSERT_DOUBLES_EQUAL(-1., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
1182   CPPUNIT_ASSERT_DOUBLES_EQUAL(-3., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
1183   CPPUNIT_ASSERT_DOUBLES_EQUAL( 5., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
1184   CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
1185   // out of range
1186   CPPUNIT_ASSERT_THROW(aSubField1->setRow(3, row), MEDEXCEPTION);
1187
1188   // setColumn
1189   double col[2] =
1190     {
1191       -7., -9.
1192     };
1193   aSubField1->setColumn(1, col);
1194   CPPUNIT_ASSERT_DOUBLES_EQUAL(-7., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
1195   CPPUNIT_ASSERT_DOUBLES_EQUAL(-3., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
1196   //#ifdef ENABLE_FORCED_FAILURES
1197   // (BUG) in MEDMEM_Array::setColumn()
1198   CPPUNIT_ASSERT_DOUBLES_EQUAL(-9., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
1199   //#endif
1200   CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
1201   // out of range
1202   CPPUNIT_ASSERT_THROW(aSubField1->setColumn(3, col), MEDEXCEPTION);
1203
1204   // setArray (Gauss)
1205   {
1206     int nbelgeoc[2] =
1207       {
1208         1, 3
1209       }; // 3 - 1 = two elements for the first (and the only) type
1210     int nbgaussgeo[2] =
1211       {
1212         -1, 1
1213       }; // one gauss point per each element
1214     MEDMEM_ArrayInterface<double,FullInterlace,Gauss>::Array * aNewArrayGauss =
1215       new MEDMEM_ArrayInterface<double,FullInterlace,Gauss>::Array
1216       (/*dim*/2, /*nbelem*/2, /*nbtypegeo*/1, /*nbelgeoc*/nbelgeoc, /*nbgaussgeo*/nbgaussgeo);
1217
1218     //#ifdef ENABLE_FAULTS
1219     aNewArrayGauss->setIJ(1, 1, -4.);
1220     aNewArrayGauss->setIJ(1, 2, -2.);
1221     aNewArrayGauss->setIJ(2, 1, -5.);
1222     aNewArrayGauss->setIJ(2, 2, -1.);
1223     //#endif
1224     //#ifdef ENABLE_FORCED_FAILURES
1225     // ? (BUG) in FullInterlaceGaussPolicy::getIndex(int i,int j)
1226     // FullInterlaceGaussPolicy::getIndex(2,2) returns 4!!!
1227     //    CPPUNIT_FAIL("? Bug in FullInterlaceGaussPolicy::getIndex(int i,int j) ?");
1228     //#endif
1229
1230     aNewArrayGauss->setIJK(1, 1, 1, -4.);
1231     aNewArrayGauss->setIJK(1, 2, 1, -2.);
1232     aNewArrayGauss->setIJK(2, 1, 1, -5.);
1233     aNewArrayGauss->setIJK(2, 2, 1, -1.);
1234
1235     aSubField1->setArray(aNewArrayGauss);
1236     // no need to delete aNewArrayGauss, because it will be deleted
1237     // in destructor or in deallocValue() method of aSubField1
1238
1239     //#ifdef ENABLE_FAULTS
1240     CPPUNIT_ASSERT_DOUBLES_EQUAL(-4., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
1241     CPPUNIT_ASSERT_DOUBLES_EQUAL(-2., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
1242     CPPUNIT_ASSERT_DOUBLES_EQUAL(-5., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
1243     CPPUNIT_ASSERT_DOUBLES_EQUAL(-1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
1244     //#endif
1245     //#ifdef ENABLE_FORCED_FAILURES
1246     // ? (BUG) in FullInterlaceGaussPolicy::getIndex(int i,int j)
1247     // Must be   : return _G[i-1]-1 + (j-1);
1248     // Instead of: return _G[i-1]-1 + (j-1)*_dim;
1249     //    CPPUNIT_FAIL("? Bug in FullInterlaceGaussPolicy::getIndex(int i,int j) ?");
1250     //#endif
1251
1252     CPPUNIT_ASSERT_DOUBLES_EQUAL(-4., aSubField1->getValueIJK(anElems1[0], 1, 1), 0.000001);
1253     CPPUNIT_ASSERT_DOUBLES_EQUAL(-2., aSubField1->getValueIJK(anElems1[0], 2, 1), 0.000001);
1254     CPPUNIT_ASSERT_DOUBLES_EQUAL(-5., aSubField1->getValueIJK(anElems1[1], 1, 1), 0.000001);
1255     CPPUNIT_ASSERT_DOUBLES_EQUAL(-1., aSubField1->getValueIJK(anElems1[1], 2, 1), 0.000001);
1256   }
1257
1258   // alloc/dealloc; compatibility of new size with support
1259   try
1260     {
1261       aSubField1->deallocValue();
1262       aSubField1->allocValue(/*NumberOfComponents*/2, /*LengthValue*/5);
1263       //#ifdef ENABLE_FAULTS
1264       // (BUG) No compatibility between Support and allocated value
1265       aSubField1->normL1();
1266       //#endif
1267       //#ifdef ENABLE_FORCED_FAILURES
1268       //    CPPUNIT_FAIL("Error: no compatibility between Support and allocated value");
1269       //#endif
1270     }
1271   catch (MEDEXCEPTION & ex)
1272     {
1273       // normal behaviour
1274     }
1275   catch (...)
1276     {
1277       CPPUNIT_FAIL("Error: no compatibility between Support and allocated value");
1278     }
1279
1280   // check that aFieldOnGroup1 is not changed after aSubField1 modifications
1281   for (int i = 1; i <= nbVals; i++)
1282     {
1283       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
1284       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
1285     }
1286
1287   // reset aFieldOnGroup2 values for simple control of operators results
1288   for (int i = 1; i <= nbVals; i++)
1289     {
1290       aFieldOnGroup2->setValueIJ(anElems[i-1], 1, i*i);
1291       aFieldOnGroup2->setValueIJ(anElems[i-1], 2, -i*i);
1292     }
1293
1294   int len = aFieldOnGroup1->getValueLength();
1295   const double * val1 = aFieldOnGroup1->getValue();
1296   const double * val2 = aFieldOnGroup2->getValue();
1297   const double * val_res;
1298
1299   // operators and add, sub, mul, div
1300
1301   // +
1302   FIELD<double> *aSum = *aFieldOnGroup1 + *aFieldOnGroup2;
1303   aSum->setName(aFieldOnGroup1->getName());
1304   aSum->setDescription(aFieldOnGroup1->getDescription());
1305   compareField_(aFieldOnGroup1, aSum, true, true);
1306   val_res = aSum->getValue();
1307   for (int i = 0; i < len; i++)
1308     {
1309       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] + val2[i], val_res[i], 0.000001);
1310     }
1311   aSum->removeReference();
1312
1313   // -
1314   FIELD<double> *aDifference = *aFieldOnGroup1 - *aFieldOnGroup2;
1315   aDifference->setName(aFieldOnGroup1->getName());
1316   aDifference->setDescription(aFieldOnGroup1->getDescription());
1317   compareField_(aFieldOnGroup1, aDifference, true, true);
1318   val_res = aDifference->getValue();
1319   for (int i = 0; i < len; i++)
1320     {
1321       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] - val2[i], val_res[i], 0.000001);
1322     }
1323   aDifference->removeReference();
1324
1325   // - (unary)
1326   FIELD<double> *aNegative = - *aFieldOnGroup1;
1327   aNegative->setName(aFieldOnGroup1->getName());
1328   aNegative->setDescription(aFieldOnGroup1->getDescription());
1329   compareField_(aFieldOnGroup1, aNegative, true, true);
1330   val_res = aNegative->getValue();
1331   for (int i = 0; i < len; i++)
1332     {
1333       CPPUNIT_ASSERT_DOUBLES_EQUAL(- val1[i], val_res[i], 0.000001);
1334     }
1335   aNegative->removeReference();
1336
1337   // *
1338   FIELD<double> *aProduct = (*aFieldOnGroup1) * (*aFieldOnGroup2);
1339   aProduct->setName(aFieldOnGroup1->getName());
1340   aProduct->setDescription(aFieldOnGroup1->getDescription());
1341   compareField_(aFieldOnGroup1, aProduct, true, true);
1342   val_res = aProduct->getValue();
1343   for (int i = 0; i < len; i++)
1344     {
1345       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] * val2[i], val_res[i], 0.000001);
1346     }
1347   aProduct->removeReference();
1348   // /
1349   FIELD<double> *aQuotient = *aFieldOnGroup1 / *aFieldOnGroup2;
1350   aQuotient->setName(aFieldOnGroup1->getName());
1351   aQuotient->setDescription(aFieldOnGroup1->getDescription());
1352   compareField_(aFieldOnGroup1, aQuotient, true, true);
1353   val_res = aQuotient->getValue();
1354   for (int i = 0; i < len; i++)
1355     {
1356       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] / val2[i], val_res[i], 0.000001);
1357     }
1358   aQuotient->removeReference();
1359
1360   double val22 = aFieldOnGroup2->getValueIJ(anElems[2], 2);
1361   aFieldOnGroup2->setValueIJ(anElems[2], 2, 0.);
1362
1363   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 / *aFieldOnGroup2, MEDEXCEPTION);
1364   //#ifdef ENABLE_FORCED_FAILURES
1365   // (BUG) is it up to user to control validity of data to avoid division on zero?
1366   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 /= *aFieldOnGroup2, MEDEXCEPTION);
1367   //#endif
1368   CPPUNIT_ASSERT_THROW(FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
1369   CPPUNIT_ASSERT_THROW(FIELD<double>::divDeep(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
1370
1371   // restore value
1372   aFieldOnGroup2->setValueIJ(anElems[2], 2, val22);
1373
1374   // restore values
1375   for (int i = 1; i <= nbVals; i++)
1376     {
1377       aFieldOnGroup1->setValueIJ(anElems[i-1], 1, 13 + i);
1378       aFieldOnGroup1->setValueIJ(anElems[i-1], 2, 13 - i);
1379     }
1380
1381   // static methods
1382   FIELD<double> * aPtr;
1383
1384   // add
1385   aPtr = FIELD<double>::add(*aFieldOnGroup1, *aFieldOnGroup2);
1386   aPtr->setName(aFieldOnGroup1->getName());
1387   aPtr->setDescription(aFieldOnGroup1->getDescription());
1388   compareField_(aFieldOnGroup1, aPtr, true, true);
1389   val_res = aPtr->getValue();
1390   for (int i = 0; i < len; i++)
1391     {
1392       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] + val2[i], val_res[i], 0.000001);
1393     }
1394   delete aPtr;
1395
1396   // sub
1397   aPtr = FIELD<double>::sub(*aFieldOnGroup1, *aFieldOnGroup2);
1398   aPtr->setName(aFieldOnGroup1->getName());
1399   aPtr->setDescription(aFieldOnGroup1->getDescription());
1400   compareField_(aFieldOnGroup1, aPtr, true, true);
1401   val_res = aPtr->getValue();
1402   for (int i = 0; i < len; i++)
1403     {
1404       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] - val2[i], val_res[i], 0.000001);
1405     }
1406   delete aPtr;
1407
1408   // mul
1409   aPtr = FIELD<double>::mul(*aFieldOnGroup1, *aFieldOnGroup2);
1410   aPtr->setName(aFieldOnGroup1->getName());
1411   aPtr->setDescription(aFieldOnGroup1->getDescription());
1412   compareField_(aFieldOnGroup1, aPtr, true, true);
1413   val_res = aPtr->getValue();
1414   for (int i = 0; i < len; i++)
1415     {
1416       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] * val2[i], val_res[i], 0.000001);
1417     }
1418   delete aPtr;
1419
1420   // div
1421   aPtr = FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup2);
1422   aPtr->setName(aFieldOnGroup1->getName());
1423   aPtr->setDescription(aFieldOnGroup1->getDescription());
1424   compareField_(aFieldOnGroup1, aPtr, true, true);
1425   val_res = aPtr->getValue();
1426   for (int i = 0; i < len; i++)
1427     {
1428       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] / val2[i], val_res[i], 0.000001);
1429     }
1430   delete aPtr;
1431
1432   // addDeep
1433   aPtr = FIELD<double>::addDeep(*aFieldOnGroup1, *aFieldOnGroup2);
1434   aPtr->setName(aFieldOnGroup1->getName());
1435   aPtr->setDescription(aFieldOnGroup1->getDescription());
1436   compareField_(aFieldOnGroup1, aPtr, true, true);
1437   val_res = aPtr->getValue();
1438   for (int i = 0; i < len; i++)
1439     {
1440       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] + val2[i], val_res[i], 0.000001);
1441     }
1442   delete aPtr;
1443
1444   // subDeep
1445   aPtr = FIELD<double>::subDeep(*aFieldOnGroup1, *aFieldOnGroup2);
1446   aPtr->setName(aFieldOnGroup1->getName());
1447   aPtr->setDescription(aFieldOnGroup1->getDescription());
1448   compareField_(aFieldOnGroup1, aPtr, true, true);
1449   val_res = aPtr->getValue();
1450   for (int i = 0; i < len; i++)
1451     {
1452       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] - val2[i], val_res[i], 0.000001);
1453     }
1454   delete aPtr;
1455
1456   // mulDeep
1457   aPtr = FIELD<double>::mulDeep(*aFieldOnGroup1, *aFieldOnGroup2);
1458   aPtr->setName(aFieldOnGroup1->getName());
1459   aPtr->setDescription(aFieldOnGroup1->getDescription());
1460   compareField_(aFieldOnGroup1, aPtr, true, true);
1461   val_res = aPtr->getValue();
1462   for (int i = 0; i < len; i++)
1463     {
1464       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] * val2[i], val_res[i], 0.000001);
1465     }
1466   delete aPtr;
1467
1468   // divDeep
1469   aPtr = FIELD<double>::divDeep(*aFieldOnGroup1, *aFieldOnGroup2);
1470   aPtr->setName(aFieldOnGroup1->getName());
1471   aPtr->setDescription(aFieldOnGroup1->getDescription());
1472   compareField_(aFieldOnGroup1, aPtr, true, true);
1473   val_res = aPtr->getValue();
1474   for (int i = 0; i < len; i++)
1475     {
1476       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] / val2[i], val_res[i], 0.000001);
1477     }
1478   delete aPtr;
1479
1480   // +=
1481   *aFieldOnGroup1 += *aFieldOnGroup2;
1482   for (int i = 1; i <= nbVals; i++)
1483     {
1484       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i + i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
1485       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i - i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
1486     }
1487
1488   // -=
1489   *aFieldOnGroup1 -= *aFieldOnGroup2;
1490   for (int i = 1; i <= nbVals; i++)
1491     {
1492       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
1493       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
1494     }
1495
1496   // *=
1497   *aFieldOnGroup1 *= *aFieldOnGroup2;
1498   for (int i = 1; i <= nbVals; i++)
1499     {
1500       CPPUNIT_ASSERT_DOUBLES_EQUAL( (13 + i)*i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
1501       CPPUNIT_ASSERT_DOUBLES_EQUAL(-(13 - i)*i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
1502     }
1503
1504   // /=
1505   *aFieldOnGroup1 /= *aFieldOnGroup2;
1506   for (int i = 1; i <= nbVals; i++)
1507     {
1508       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
1509       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
1510     }
1511
1512   // check case of different operands: support
1513   MESH * aMeshOneMore = MEDMEMTest_createTestMesh();
1514   const GROUP* aGroupOneMore = aMeshOneMore->getGroup(MED_EN::MED_FACE, 1);
1515   FIELD<double> * aFieldOnGroup3 =
1516     createFieldOnGroup<double>(aMeshOneMore, aGroupOneMore, "Test_Diff_Mesh", "test");
1517   for (int i = 1; i <= nbVals; i++)
1518     {
1519       aFieldOnGroup3->setValueIJ(anElems[i-1], 1, 2*i);
1520       aFieldOnGroup3->setValueIJ(anElems[i-1], 2, 3*i);
1521     }
1522   const double * val3 = aFieldOnGroup3->getValue();
1523
1524   //CPPUNIT_ASSERT_NO_THROW();
1525   try
1526     {
1527       aPtr = FIELD<double>::addDeep(*aFieldOnGroup1, *aFieldOnGroup3);
1528       aPtr->setName(aFieldOnGroup1->getName());
1529       aPtr->setDescription(aFieldOnGroup1->getDescription());
1530       compareField_(aFieldOnGroup1, aPtr, true, true);
1531       val_res = aPtr->getValue();
1532       for (int i = 0; i < len; i++)
1533         {
1534           CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] + val3[i], val_res[i], 0.000001);
1535         }
1536       delete aPtr;
1537
1538       aPtr = FIELD<double>::subDeep(*aFieldOnGroup1, *aFieldOnGroup3);
1539       delete aPtr;
1540       aPtr = FIELD<double>::mulDeep(*aFieldOnGroup1, *aFieldOnGroup3);
1541       delete aPtr;
1542       aPtr = FIELD<double>::divDeep(*aFieldOnGroup1, *aFieldOnGroup3);
1543       delete aPtr;
1544     }
1545   catch (MEDEXCEPTION & ex)
1546     {
1547       CPPUNIT_FAIL(ex.what());
1548     }
1549   catch (...)
1550     {
1551       CPPUNIT_FAIL("Unknown exception in FIELD::xxxDeep()");
1552     }
1553
1554   CPPUNIT_ASSERT_THROW(FIELD<double>::add(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
1555   CPPUNIT_ASSERT_THROW(FIELD<double>::sub(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
1556   CPPUNIT_ASSERT_THROW(FIELD<double>::mul(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
1557   CPPUNIT_ASSERT_THROW(FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
1558
1559   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 + *aFieldOnGroup3, MEDEXCEPTION);
1560   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 - *aFieldOnGroup3, MEDEXCEPTION);
1561   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 * *aFieldOnGroup3, MEDEXCEPTION);
1562   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 / *aFieldOnGroup3, MEDEXCEPTION);
1563
1564   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 += *aFieldOnGroup3, MEDEXCEPTION);
1565   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 -= *aFieldOnGroup3, MEDEXCEPTION);
1566   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 *= *aFieldOnGroup3, MEDEXCEPTION);
1567   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 /= *aFieldOnGroup3, MEDEXCEPTION);
1568
1569   // check case of different operands: MEDComponentsUnits
1570   aFieldOnGroup1->setMEDComponentUnit(1, "unit3");
1571
1572   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 + *aFieldOnGroup2, MEDEXCEPTION);
1573   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 - *aFieldOnGroup2, MEDEXCEPTION);
1574   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 += *aFieldOnGroup2, MEDEXCEPTION);
1575   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 -= *aFieldOnGroup2, MEDEXCEPTION);
1576   CPPUNIT_ASSERT_THROW(FIELD<double>::add(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
1577   CPPUNIT_ASSERT_THROW(FIELD<double>::sub(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
1578   CPPUNIT_ASSERT_THROW(FIELD<double>::addDeep(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
1579   CPPUNIT_ASSERT_THROW(FIELD<double>::subDeep(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
1580
1581   //CPPUNIT_ASSERT_NO_THROW();
1582   try
1583     {
1584       aPtr = FIELD<double>::mul(*aFieldOnGroup1, *aFieldOnGroup2);
1585       delete aPtr;
1586       aPtr = FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup2);
1587       delete aPtr;
1588       aPtr = FIELD<double>::mulDeep(*aFieldOnGroup1, *aFieldOnGroup2);
1589       delete aPtr;
1590       aPtr = FIELD<double>::divDeep(*aFieldOnGroup1, *aFieldOnGroup2);
1591       delete aPtr;
1592
1593       *aFieldOnGroup1 *= *aFieldOnGroup2;
1594       *aFieldOnGroup1 /= *aFieldOnGroup2;
1595
1596       FIELD<double> *aPr = *aFieldOnGroup1 * *aFieldOnGroup2;
1597       FIELD<double> *aQu = *aFieldOnGroup1 / *aFieldOnGroup2;
1598       aPr->removeReference();
1599       aQu->removeReference();
1600     }
1601   catch (MEDEXCEPTION & ex)
1602     {
1603       CPPUNIT_FAIL(ex.what());
1604     }
1605   catch (...)
1606     {
1607       CPPUNIT_FAIL("Unknown exception");
1608     }
1609
1610   // restore MED units
1611   aFieldOnGroup1->setMEDComponentUnit(1, "unit1");
1612
1613   // check case of different operands: numberOfComponents
1614   //#ifdef ENABLE_FAULTS
1615   // (BUG) Cannot allocate value of higher dimension because of _componentsTypes reinitialization
1616   // Must be MEDEXCEPTION instead. And on attempt to change nb.components must be the same behaviour.
1617   aFieldOnGroup1->deallocValue();
1618   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->allocValue(/*dim*/5), MEDEXCEPTION);
1619   //#endif
1620
1621   aSubSupport1->removeReference();
1622   delete [] anElems1;
1623
1624   delete aScalarProduct;
1625   delete aSubField1;
1626   delete anAreaField;
1627   delete barycenter;
1628   delete aFieldOnGroup1;
1629   delete aFieldOnGroup2;
1630   delete aFieldOnGroup3;
1631
1632   delete aMesh;
1633   delete aMeshOneMore;
1634
1635   /////////////////////
1636   // TEST 5: Drivers //
1637   /////////////////////
1638   testDrivers();
1639 }
1640
1641 int main (int argc, char** argv)
1642 {
1643   MEDMEMTest_testField();
1644 }