src/MEDMEMCppTest/MEDMEMTest_Field.cxx

   1 // Copyright (C) 2007-2012  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_MedMeshDriver.hxx"
  25 #include "MEDMEM_Grid.hxx"
  26 #include "MEDMEM_Group.hxx"
  27 #include "MEDMEM_Mesh.hxx"
  28 #include "MEDMEM_Meshing.hxx"
  29 #include "MEDMEM_Support.hxx"
  30 #include "MEDMEM_VtkMeshDriver.hxx"
  31
  32
  33 #include <cppunit/TestAssert.h>
  34 #include <sstream>
  35 #include <cmath>
  36
  37 // use this define to enable lines, execution of which leads to Segmentation Fault
  38 //#define ENABLE_FAULTS
  39
  40 // use this define to enable CPPUNIT asserts and fails, showing bugs
  41 //#define ENABLE_FORCED_FAILURES
  42
  43 using namespace std;
  44 using namespace MEDMEM;
  45
  46 // #14,15: MEDMEMTest_Field.cxx
  47 // Check methods from MEDMEM_Field.hxx, MEDMEM_FieldConvert.hxx
  48
  49 /*!
  50  *  Check methods (48), defined in MEDMEM_Field.hxx:
  51  *  class FIELD_
  52 {
  53  *   (+)     FIELD_();
  54  *   (+)     FIELD_(const SUPPORT * Support, const int NumberOfComponents);
  55  *   (+)     FIELD_(const FIELD_ &m);
  56  *   (+)     virtual ~FIELD_();
  57  *   (+)     FIELD_& operator=(const FIELD_ &m);
  58  *
  59  *   (-)     virtual  void     rmDriver(int index=0);
  60  *   (-)     virtual   int     addDriver(driverTypes driverType,
  61  *                             const string & fileName="Default File Name.med",
  62  *                              const string & driverFieldName="Default Field Nam",
  63  *                              MED_EN::med_mode_acces access=MED_EN::MED_REMP);
  64  *   (-)     virtual  int      addDriver(GENDRIVER & driver);
  65  *
  66  *   (-)     virtual  void     read (const GENDRIVER &);
  67  *   (-)     virtual  void     read(int index=0);
  68  *   (-)     virtual  void     openAppend(void);
  69  *   (-)     virtual  void     write(const GENDRIVER &);
  70  *   (-)     virtual  void     write(int index=0, const string & driverName="");
  71  *   (-)     virtual  void     writeAppend(const GENDRIVER &);
  72  *   (-)     virtual  void     writeAppend(int index=0, const string & driverName="");
  73  *
  74  *   (+)     inline void     setName(const string Name);
  75  *   (+)     inline string   getName() const;
  76  *   (+)     inline void     setDescription(const string Description);
  77  *   (+)     inline string   getDescription() const;
  78  *   (+)     inline const SUPPORT * getSupport() const;
  79  *   (+)     inline void     setSupport(const SUPPORT * support);
  80  *   (+)     inline void     setNumberOfComponents(const int NumberOfComponents);
  81  *   (+)     inline int      getNumberOfComponents() const;
  82  *   (+)     inline void     setNumberOfValues(const int NumberOfValues);
  83  *   (+)     inline int      getNumberOfValues() const;
  84  *   (+)     inline void     setComponentsNames(const string * ComponentsNames);
  85  *   (+)     inline void     setComponentName(int i, const string ComponentName);
  86  *   (+)     inline const string * getComponentsNames() const;
  87  *   (+)     inline string   getComponentName(int i) const;
  88  *   (+)     inline void     setComponentsDescriptions(const string * ComponentsDescriptions);
  89  *   (+)     inline void     setComponentDescription(int i, const string ComponentDescription);
  90  *   (+)     inline const string * getComponentsDescriptions() const;
  91  *   (+)     inline string   getComponentDescription(int i) const;
  92  *   (+)     inline void     setComponentsUnits(const UNIT * ComponentsUnits);
  93  *   (+)     inline const UNIT *   getComponentsUnits() const;
  94  *   (+)     inline const UNIT *   getComponentUnit(int i) const;
  95  *   (+)     inline void     setMEDComponentsUnits(const string * MEDComponentsUnits);
  96  *   (+)     inline void     setMEDComponentUnit(int i, const string MEDComponentUnit);
  97  *   (+)     inline const string * getMEDComponentsUnits() const;
  98  *   (+)     inline string   getMEDComponentUnit(int i) const;
  99  *
 100  *   (+)     inline void     setIterationNumber(int IterationNumber);
 101  *   (+)     inline int      getIterationNumber() const;
 102  *   (+)     inline void     setTime(double Time);
 103  *   (+)     inline double   getTime() const;
 104  *   (+)     inline void     setOrderNumber(int OrderNumber);
 105  *   (+)     inline int      getOrderNumber() const;
 106  *
 107  *   (+)     inline MED_EN::med_type_champ getValueType () const;
 108  *   (+)     inline MED_EN::medModeSwitch  getInterlacingType() const;
 109  *   (-)     virtual inline bool getGaussPresence() const throw (MEDEXCEPTION);
 110  *
 111 }
 112  *
 113  *  template <class T, class INTERLACING_TAG> class FIELD : public FIELD_
 114 {
 115  *   (+)     FIELD();
 116  *   (+)     FIELD(const FIELD &m);
 117  *   (+)     FIELD(const SUPPORT * Support, const int NumberOfComponents) throw (MEDEXCEPTION);
 118  *   (+)     FIELD(driverTypes driverType,
 119  *                 const string & fileName, const string & fieldDriverName,
 120  *                 const int iterationNumber=-1, const int orderNumber=-1) throw (MEDEXCEPTION);
 121  *   (+) FIELD(const SUPPORT * Support, driverTypes driverType,
 122  *                 const string & fileName="", const string & fieldName="",
 123  *                 const int iterationNumber = -1, const int orderNumber = -1) throw (MEDEXCEPTION);
 124  *   (+)     ~FIELD();
 125  *   (+)     FIELD & operator=(const FIELD &m);
 126  *
 127  *   (+) const FIELD operator+(const FIELD& m) const;
 128  *   (+) const FIELD operator-(const FIELD& m) const;
 129  *   (+) const FIELD operator*(const FIELD& m) const;
 130  *   (+) const FIELD operator/(const FIELD& m) const;
 131  *   (+) const FIELD operator-() const;
 132  *   (+) FIELD& operator+=(const FIELD& m);
 133  *   (+) FIELD& operator-=(const FIELD& m);
 134  *   (+) FIELD& operator*=(const FIELD& m);
 135  *   (+) FIELD& operator/=(const FIELD& m);
 136  *
 137  *   (+) static FIELD* add(const FIELD& m, const FIELD& n);
 138  *   (+) static FIELD* addDeep(const FIELD& m, const FIELD& n);
 139  *   (+) static FIELD* sub(const FIELD& m, const FIELD& n);
 140  *   (+) static FIELD* subDeep(const FIELD& m, const FIELD& n);
 141  *   (+) static FIELD* mul(const FIELD& m, const FIELD& n);
 142  *   (+) static FIELD* mulDeep(const FIELD& m, const FIELD& n);
 143  *   (+) static FIELD* div(const FIELD& m, const FIELD& n);
 144  *   (+) static FIELD* divDeep(const FIELD& m, const FIELD& n);
 145  *
 146  *   (+)     double normMax() const throw (MEDEXCEPTION);
 147  *   (+)     double norm2() const throw (MEDEXCEPTION);
 148  *
 149  *   (+)     void   applyLin(T a, T b);
 150  *   (+)     template <T T_function(T)> void applyFunc();
 151  *   (+)     void applyPow(T scalar);
 152  *
 153  *   (+)     static FIELD* scalarProduct(const FIELD& m, const FIELD& n, bool deepCheck=false);
 154  *
 155  *   (+)     double normL2(int component, const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 156  *   (+)     double normL2(const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 157  *   (+)     double normL1(int component, const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 158  *   (+)     double normL1(const FIELD<double,FullInterlace> * p_field_volume=NULL) const;
 159  *
 160  *   (+)     FIELD* extract(const SUPPORT *subSupport) const throw (MEDEXCEPTION);
 161  *
 162  *   (EMPTY COMMENT, EMPTY IMPLEMENTATION!!!) void init ();
 163  *
 164  *   (+)     void rmDriver(int index=0);
 165  *   (+)     int  addDriver(driverTypes driverType,
 166  *                          const string & fileName="Default File Name.med",
 167  *                          const string & driverFieldName="Default Field Name",
 168  *                          MED_EN::med_mode_acces access=MED_EN::MED_REMP);
 169  *   (+)     int  addDriver(GENDRIVER & driver);
 170  *
 171  *   (+)     void allocValue(const int NumberOfComponents);
 172  *   (+)     void allocValue(const int NumberOfComponents, const int LengthValue);
 173  *   (+)     void deallocValue();
 174  *
 175  *   (+)     inline void read(int index=0);
 176  *   (+)     inline void read(const GENDRIVER & genDriver);
 177  *   (+)     inline void read(driverTypes driverType, const std::string& filename);
 178  *   (+)     inline void write(int index=0);
 179  *   (+)     inline void write(const GENDRIVER &,
 180  *                             MED_EN::med_mode_acces medMode=MED_EN::RDWR);
 181  *   (+)     inline void write(driverTypes driverType, const std::string& filename,
 182  *                             MED_EN::med_mode_acces medMode=MED_EN::RDWR);
 183  *   (+)     inline void writeAppend(int index=0, const string & driverName = "");
 184  *   (+) inline void writeAppend(const GENDRIVER &);
 185  *
 186  *   (+)     inline MEDMEM_Array_  * getArray()        const throw (MEDEXCEPTION);
 187  *   (+)     inline ArrayGauss     * getArrayGauss()   const throw (MEDEXCEPTION);
 188  *   (+)     inline ArrayNoGauss   * getArrayNoGauss() const throw (MEDEXCEPTION);
 189  *   (+)     inline bool             getGaussPresence() const throw (MEDEXCEPTION);
 190  *
 191  *   (+)     inline int      getValueLength() const throw (MEDEXCEPTION);
 192  *   (+)     inline const T* getValue()       const throw (MEDEXCEPTION);
 193  *   (+)     inline const T* getRow(int i)    const throw (MEDEXCEPTION);
 194  *   (+)     inline const T* getColumn(int j) const throw (MEDEXCEPTION);
 195  *   (+)     inline T        getValueIJ(int i,int j) const throw (MEDEXCEPTION);
 196  *   (+)     inline T        getValueIJK(int i,int j,int k) const throw (MEDEXCEPTION);
 197  *   (+)     bool            getValueOnElement(int eltIdInSup,T* retValues) const throw (MEDEXCEPTION);
 198  *
 199  *   (+)     const int getNumberOfGeometricTypes() const throw (MEDEXCEPTION);
 200  *
 201  *   (+)     const GAUSS_LOCALIZATION<INTERLACING_TAG> & getGaussLocalization
 202  *                          (MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
 203  *   (+)     const GAUSS_LOCALIZATION<INTERLACING_TAG> * getGaussLocalizationPtr
 204  *                          (MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
 205  *   (+)     void setGaussLocalization(MED_EN::medGeometryElement geomElement,
 206  *                                     const GAUSS_LOCALIZATION<INTERLACING_TAG> & gaussloc);
 207  *   (+)     const int * getNumberOfGaussPoints() const throw (MEDEXCEPTION);
 208  *   (+)     const int   getNumberOfGaussPoints
 209  *                          (MED_EN::medGeometryElement geomElement) const throw (MEDEXCEPTION);
 210  *   (+)     const int   getNbGaussI(int i) const throw (MEDEXCEPTION);
 211  *
 212  *   (+)     const int * getNumberOfElements() const throw (MEDEXCEPTION);
 213  *   (+)     const MED_EN::medGeometryElement * getGeometricTypes() const throw (MEDEXCEPTION);
 214  *   (+)     bool        isOnAllElements() const throw (MEDEXCEPTION);
 215  *
 216  *   (+)     inline void setArray(MEDMEM_Array_ *value) throw (MEDEXCEPTION);
 217  *   (+)     FIELD<double, FullInterlace>* getGaussPointsCoordinates() const throw (MEDEXCEPTION);
 218  *   (+)     inline void setValue(T* value) throw (MEDEXCEPTION);
 219  *   (+)     inline void setRow(int i, T* value) throw (MEDEXCEPTION);
 220  *   (+)     inline void setColumn(int i, T* value) throw (MEDEXCEPTION);
 221  *   (+)     inline void setValueIJ(int i, int j, T value) throw (MEDEXCEPTION);
 222  *
 223  *   (NOT IMPLEMENTED!!!) void getVolume() const throw (MEDEXCEPTION);
 224  *   (NOT IMPLEMENTED!!!) void getArea() const throw (MEDEXCEPTION);
 225  *   (NOT IMPLEMENTED!!!) void getLength() const throw (MEDEXCEPTION);
 226  *   (NOT IMPLEMENTED!!!) void getNormal() const throw (MEDEXCEPTION);
 227  *   (NOT IMPLEMENTED!!!) void getBarycenter() const throw (MEDEXCEPTION);
 228  *
 229  *   (+)     void fillFromAnalytic(myFuncType f) throw (MEDEXCEPTION);
 230  *
 231 }
 232  *
 233  *  Use code of test_operation_fieldint.cxx
 234  *              test_operation_fielddouble.cxx
 235  *              test_copie_field_.cxx
 236  *              test_copie_fieldT.cxx
 237  */
 238  static void compareField_(const FIELD_ * theField_1, const FIELD_ * theField_2, bool isFIELD, bool isValue)
 239 {
 240   // name, description, support
 241   CPPUNIT_ASSERT_EQUAL(theField_1->getName(), theField_2->getName());
 242   CPPUNIT_ASSERT_EQUAL(theField_1->getDescription(), theField_2->getDescription());
 243   CPPUNIT_ASSERT_EQUAL(theField_1->getSupport(), theField_2->getSupport());
 244
 245   // components information
 246   int aNbComps = theField_1->getNumberOfComponents();
 247   CPPUNIT_ASSERT_EQUAL(aNbComps, theField_2->getNumberOfComponents());
 248
 249   for (int i = 1; i <= aNbComps; i++)
 250     {
 251       CPPUNIT_ASSERT_EQUAL(theField_1->getComponentName(i), theField_2->getComponentName(i));
 252       CPPUNIT_ASSERT_EQUAL(theField_1->getComponentDescription(i), theField_2->getComponentDescription(i));
 253       CPPUNIT_ASSERT_EQUAL(theField_1->getMEDComponentUnit(i), theField_2->getMEDComponentUnit(i));
 254     }
 255
 256   // iteration information
 257   CPPUNIT_ASSERT_EQUAL(theField_1->getIterationNumber(), theField_2->getIterationNumber());
 258   CPPUNIT_ASSERT_EQUAL(theField_1->getOrderNumber(), theField_2->getOrderNumber());
 259   CPPUNIT_ASSERT_DOUBLES_EQUAL(theField_1->getTime(), theField_2->getTime(), 0.0000001);
 260
 261   // Value
 262   int nbOfValues = theField_1->getNumberOfValues();
 263   CPPUNIT_ASSERT_EQUAL(nbOfValues, theField_2->getNumberOfValues());
 264
 265   if (isFIELD)
 266     {
 267       // Value type and Interlacing type
 268       CPPUNIT_ASSERT_EQUAL(theField_1->getValueType(), theField_2->getValueType());
 269       CPPUNIT_ASSERT_EQUAL(theField_1->getInterlacingType(), theField_2->getInterlacingType());
 270
 271       // Gauss Presence
 272       if (isValue)
 273         {
 274           CPPUNIT_ASSERT_EQUAL(theField_1->getGaussPresence(), theField_2->getGaussPresence());
 275         }
 276       else
 277         {
 278           CPPUNIT_ASSERT_THROW(theField_1->getGaussPresence(), MEDEXCEPTION);
 279           CPPUNIT_ASSERT_THROW(theField_2->getGaussPresence(), MEDEXCEPTION);
 280         }
 281     }
 282   else
 283     {
 284       CPPUNIT_ASSERT_THROW(theField_1->getGaussPresence(), MEDEXCEPTION);
 285       CPPUNIT_ASSERT_THROW(theField_2->getGaussPresence(), MEDEXCEPTION);
 286     }
 287 }
 288
 289 static void checkField_(FIELD_ * theField_, const SUPPORT * theSupport,
 290                         MED_EN::med_type_champ theValueType,
 291                         MED_EN::medModeSwitch theInterlace)
 292 {
 293   // name
 294   const string aFieldName = "a_name_of_a_field";
 295   theField_->setName(aFieldName);
 296   CPPUNIT_ASSERT_EQUAL(aFieldName, theField_->getName());
 297
 298   // description
 299   const string aFieldDescr = "a_description_of_a_field";
 300   theField_->setDescription(aFieldDescr);
 301   CPPUNIT_ASSERT_EQUAL(aFieldDescr, theField_->getDescription());
 302
 303   // support
 304   theField_->setSupport(theSupport);
 305   CPPUNIT_ASSERT(theField_->getSupport() == theSupport);
 306
 307   // components information
 308   int aNbComps = 3;
 309
 310   string aCompsNames[3] =
 311     {
 312       "Vx", "Vy", "Vz"
 313     };
 314   string aCompsDescs[3] =
 315     {
 316       "vitesse selon x", "vitesse selon y", "vitesse selon z"
 317     };
 318   string aCompsUnits[3] =
 319     {
 320       "m.s-1", "m.s-1", "m.s-1"
 321     };
 322
 323   theField_->setNumberOfComponents(aNbComps);
 324   CPPUNIT_ASSERT_EQUAL(aNbComps, theField_->getNumberOfComponents());
 325
 326   theField_->setComponentsNames(aCompsNames);
 327
 328   //#ifdef ENABLE_FAULTS
 329   try
 330     {
 331       theField_->setNumberOfComponents(7);
 332       // Segmentation fault here because array of components names is not resized
 333       for (int i = 1; i <= 7; i++)
 334         {
 335           theField_->setComponentName(i, "AnyComponent");
 336         }
 337     }
 338   catch (MEDEXCEPTION& ex)
 339     {
 340       // Ok, it is good to have MEDEXCEPTION here
 341     }
 342   catch (...)
 343     {
 344       CPPUNIT_FAIL("Unknown exception cought");
 345     }
 346   // restore components names
 347   theField_->setNumberOfComponents(aNbComps);
 348   theField_->setComponentsNames(aCompsNames);
 349
 350   theField_->setComponentsDescriptions(aCompsDescs);
 351   theField_->setMEDComponentsUnits(aCompsUnits);
 352
 353   const string * aCompsNamesBack = theField_->getComponentsNames();
 354   const string * aCompsDescsBack = theField_->getComponentsDescriptions();
 355   const string * aCompsUnitsBack = theField_->getMEDComponentsUnits();
 356   for (int i = 1; i <= aNbComps; i++)
 357     {
 358       CPPUNIT_ASSERT_EQUAL(aCompsNamesBack[i-1], theField_->getComponentName(i));
 359       CPPUNIT_ASSERT_EQUAL(aCompsNamesBack[i-1], aCompsNames[i-1]);
 360
 361       CPPUNIT_ASSERT_EQUAL(aCompsDescsBack[i-1], theField_->getComponentDescription(i));
 362       CPPUNIT_ASSERT_EQUAL(aCompsDescsBack[i-1], aCompsDescs[i-1]);
 363
 364       CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack[i-1], theField_->getMEDComponentUnit(i));
 365       CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack[i-1], aCompsUnits[i-1]);
 366     }
 367
 368   const string aCompName2 ("Name of second component");
 369   const string aCompDesc2 ("Description of second component");
 370   const string aCompUnit2 ("Unit of second MED component");
 371
 372   theField_->setComponentName(2, aCompName2);
 373   theField_->setComponentDescription(2, aCompDesc2);
 374   theField_->setMEDComponentUnit(2, aCompUnit2);
 375
 376   const string * aCompsNamesBack2 = theField_->getComponentsNames();
 377   const string * aCompsDescsBack2 = theField_->getComponentsDescriptions();
 378   const string * aCompsUnitsBack2 = theField_->getMEDComponentsUnits();
 379
 380   CPPUNIT_ASSERT_EQUAL(aCompsNamesBack2[1], theField_->getComponentName(2));
 381   CPPUNIT_ASSERT_EQUAL(aCompsNamesBack2[1], aCompName2);
 382
 383   CPPUNIT_ASSERT_EQUAL(aCompsDescsBack2[1], theField_->getComponentDescription(2));
 384   CPPUNIT_ASSERT_EQUAL(aCompsDescsBack2[1], aCompDesc2);
 385
 386   CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack2[1], theField_->getMEDComponentUnit(2));
 387   CPPUNIT_ASSERT_EQUAL(aCompsUnitsBack2[1], aCompUnit2);
 388
 389   CPPUNIT_ASSERT_THROW(theField_->setComponentName(0, "str"), MEDEXCEPTION);
 390   CPPUNIT_ASSERT_THROW(theField_->setComponentName(aNbComps + 1, "str"), MEDEXCEPTION);
 391   CPPUNIT_ASSERT_THROW(theField_->setComponentDescription(0, "str"), MEDEXCEPTION);
 392   CPPUNIT_ASSERT_THROW(theField_->setComponentDescription(aNbComps + 1, "str"), MEDEXCEPTION);
 393   CPPUNIT_ASSERT_THROW(theField_->setMEDComponentUnit(0, "str"), MEDEXCEPTION);
 394   CPPUNIT_ASSERT_THROW(theField_->setMEDComponentUnit(aNbComps + 1, "str"), MEDEXCEPTION);
 395
 396   // iteration information
 397   int anIterNumber = 10; // set value to MED_NOPDT if undefined (default)
 398   theField_->setIterationNumber(anIterNumber);
 399   CPPUNIT_ASSERT_EQUAL(anIterNumber, theField_->getIterationNumber());
 400
 401   int anOrderNumber = 1; // set value to MED_NONOR if undefined (default)
 402   theField_->setOrderNumber(anOrderNumber);
 403   CPPUNIT_ASSERT_EQUAL(anOrderNumber, theField_->getOrderNumber());
 404
 405   double aTime = 3.435678; // in second
 406   theField_->setTime(aTime);
 407   CPPUNIT_ASSERT_DOUBLES_EQUAL(aTime, theField_->getTime(), 0.0000001);
 408
 409   // Value
 410   int nbOfValues = 10;
 411   // dangerous method, because it does not reallocate values array
 412   theField_->setNumberOfValues(nbOfValues);
 413   CPPUNIT_ASSERT_EQUAL(nbOfValues, theField_->getNumberOfValues());
 414
 415   // Value type and Interlacing type
 416   CPPUNIT_ASSERT_EQUAL(theValueType, theField_->getValueType());
 417   CPPUNIT_ASSERT_EQUAL(theInterlace, theField_->getInterlacingType());
 418 }
 419
 420 template<class T, class INTERLACING_TAG>
 421 void compareField(const FIELD<T, INTERLACING_TAG> * theField_1,
 422                   const FIELD<T, INTERLACING_TAG> * theField_2, bool isValue)
 423 {
 424   // compare FIELD_ part
 425   compareField_(theField_1, theField_2, /*isFIELD = */true, isValue);
 426
 427   // compare FIELD part
 428   // TO DO
 429 }
 430
 431 template<class T, class INTERLACING_TAG>
 432 void checkField (FIELD<T, INTERLACING_TAG> * theField, const SUPPORT * theSupport)
 433 {
 434   // check FIELD_ part
 435   MED_EN::med_type_champ aValueType = SET_VALUE_TYPE<T>::_valueType;
 436   MED_EN::medModeSwitch  anInterlace = SET_INTERLACING_TYPE<INTERLACING_TAG>::_interlacingType;
 437   checkField_(theField, theSupport, aValueType, anInterlace);
 438
 439   // check FIELD part
 440
 441   // filling by support charackteristics (NOT IMPLEMENTED METHODS!!!):
 442   // value type must be MED_REEL64 (i.e. a FIELD<double>) for these methods,
 443   // nb. of components must be equal 1 (for Volume, Area, Length) or
 444   // space dimension (for Normal, Barycenter, )
 445   {
 446     const GMESH* aMesh = theSupport->getMesh();
 447     int spaceDim = 3;
 448     if (aMesh) spaceDim = aMesh->getSpaceDimension();
 449     theField->deallocValue();
 450     theField->allocValue(/*NumberOfComponents = */spaceDim + 1);
 451
 452     theField->deallocValue();
 453     theField->allocValue(/*NumberOfComponents = */1);
 454
 455     if (aMesh)
 456       {
 457         theField->deallocValue();
 458         theField->allocValue(/*NumberOfComponents = */spaceDim);
 459       }
 460
 461     if (aMesh)
 462       {
 463         theField->deallocValue();
 464         theField->allocValue(/*NumberOfComponents = */spaceDim);
 465         //  0020142: [CEA 315] Unused function in MEDMEM::FIELD
 466         // getVolume() etc. does nothing
 467         //       if (aValueType == MED_EN::MED_REEL64)
 468         {
 469           //         CPPUNIT_ASSERT_NO_THROW(theField->getNormal());
 470           //         CPPUNIT_ASSERT_NO_THROW(theField->getBarycenter());
 471           //
 472         }
 473         //       else
 474         {
 475           //         CPPUNIT_ASSERT_THROW(theField->getNormal(), MEDEXCEPTION);
 476           //         CPPUNIT_ASSERT_THROW(theField->getBarycenter(), MEDEXCEPTION);
 477           //
 478         }
 479       }
 480   }
 481
 482   // values
 483   theField->deallocValue();
 484   theField->allocValue(/*NumberOfComponents = */2);
 485   int nbElemSupport = theSupport->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
 486   CPPUNIT_ASSERT_EQUAL(nbElemSupport, theField->getNumberOfValues());
 487
 488   theField->deallocValue();
 489   CPPUNIT_ASSERT_THROW(theField->getGaussPresence(), MEDEXCEPTION);
 490
 491   // copy constructor
 492   FIELD<T, INTERLACING_TAG> *aField_copy1= new FIELD<T, INTERLACING_TAG>(*theField);
 493   compareField(theField, aField_copy1, /*isValue = */false);
 494   aField_copy1->removeReference();
 495
 496   // operator=
 497   FIELD<T, INTERLACING_TAG> *aField_copy2=new FIELD<T, INTERLACING_TAG>();
 498   *aField_copy2 = *theField;
 499   compareField(theField, aField_copy2, /*isValue = */false);
 500   aField_copy2->removeReference();
 501 }
 502
 503 template<class T>
 504 FIELD<T> * createFieldOnGroup(MESH* theMesh, const GROUP* theGroup,
 505                               const string theName, const string theDescr)
 506 {
 507   FIELD<T> * aFieldOnGroup = new FIELD<T>(theGroup, /*NumberOfComponents = */2);
 508
 509   aFieldOnGroup->setName(theName);
 510   aFieldOnGroup->setDescription(theDescr);
 511
 512   string aCompsNames[2] =
 513     {
 514       "Pos", "Neg"
 515     };
 516   string aCompsDescs[2] =
 517     {
 518       "+", "-"
 519     };
 520   string aCompsUnits[2] =
 521     {
 522       "unit1", "unit2"
 523     };
 524
 525   aFieldOnGroup->setComponentsNames(aCompsNames);
 526   aFieldOnGroup->setComponentsDescriptions(aCompsDescs);
 527   aFieldOnGroup->setMEDComponentsUnits(aCompsUnits);
 528
 529   return aFieldOnGroup;
 530 }
 531
 532 double plus13 (double val);
 533 double plus13 (double val)
 534 {
 535   return val + 13;
 536 }
 537
 538 // function to calculate field values from coordinates of an element
 539 // typedef void (*myFuncType)(const double * temp, T* output);
 540 // size of temp array = space dim = 3
 541 // size of output array = nb. comps = 2
 542 static void proj2d (const double * temp, double* output)
 543 {
 544   // dimetric projection with coefficients:
 545   // 1.0 along Oy and Oz, 0.5 along Ox
 546   //
 547   //    ^ z (y_)
 548   //    |
 549   //    |
 550   //    .----> y (x_)
 551   //   /
 552   //  L x
 553   //
 554   // x_ = y - x * sqrt(2.) / 4.
 555   // y_ = z - x * sqrt(2.) / 4.
 556
 557   double dx = temp[0] * std::sqrt(2.) / 4.;
 558   output[0] = temp[1] - dx;
 559   output[1] = temp[2] - dx;
 560 }
 561
 562 static void testDrivers()
 563 {
 564   string filename_rd                  = getResourceFile("pointe.med");
 565   string filename_wr                  = makeTmpFile("myMedFieldfile.med", filename_rd);
 566   string filename_support_wr          = makeTmpFile("myMedSupportFiledfile.med");
 567   string filename22_rd                = getResourceFile("pointe.med");
 568   string filenamevtk_wr               = makeTmpFile("myMedFieldfile22.vtk");
 569
 570   string fieldname_celldouble_rd      = "fieldcelldoublescalar";
 571   string fieldname_celldouble_wr      = fieldname_celldouble_rd + "_cpy";
 572   string fieldname_nodeint_rd         = "fieldnodeint";
 573   string fieldname_nodeint_wr         = fieldname_nodeint_rd + "_cpy";
 574   string fieldname_nodeint_wr1        = fieldname_nodeint_rd + "_cpy1";
 575   string meshname                     = "maa1";
 576
 577   // To remove tmp files from disk
 578   MEDMEMTest_TmpFilesRemover aRemover;
 579   aRemover.Register(filename_wr);
 580   aRemover.Register(filenamevtk_wr);
 581   aRemover.Register(filename_support_wr);
 582
 583   FIELD<int> *aInvalidField=new FIELD<int>();
 584   //must throw becase only VTK_DRIVER or MED_DRIVER may be specified as driverType for FIELD
 585   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(NO_DRIVER, filename_rd, fieldname_nodeint_rd)),
 586                        MEDEXCEPTION);
 587   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(GIBI_DRIVER, filename_rd, fieldname_nodeint_rd)),
 588                        MEDEXCEPTION);
 589   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(PORFLOW_DRIVER, filename_rd, fieldname_nodeint_rd)),
 590                        MEDEXCEPTION);
 591   CPPUNIT_ASSERT_THROW(*aInvalidField = *(new FIELD<int>(ASCII_DRIVER, filename_rd, fieldname_nodeint_rd)),
 592                        MEDEXCEPTION);
 593   aInvalidField->removeReference();
 594   //////////////////
 595   //TestRead Part//
 596   //////////////////
 597   FIELD<double> *aField_1 = NULL;
 598   CPPUNIT_ASSERT_NO_THROW(aField_1 = new FIELD<double>(MED_DRIVER, filename_rd,
 599                                                        fieldname_celldouble_rd));
 600
 601   //Test read(int index) method
 602   int IdDriver_rd = aField_1->addDriver(MED_DRIVER,filename_rd,fieldname_celldouble_rd,MED_EN::RDONLY);
 603   // TODO: throw if read for the second time
 604   // (BUG) Cannot open file, but file exist
 605   // EAP: no more pb with opening the file for the second time with "weaker" mode,
 606   // but why to re-read the field?
 607   //CPPUNIT_ASSERT_THROW(aField_1->read(IdDriver_rd),MEDEXCEPTION);
 608   CPPUNIT_ASSERT_NO_THROW(aField_1->read(IdDriver_rd));
 609
 610   //Test read(GENDRIVER & genDriver) method
 611   FIELD<int> *aField_2 = new FIELD<int>();
 612   aField_2->setName(fieldname_nodeint_rd);
 613   {
 614     MED_FIELD_RDONLY_DRIVER<int> aMedRdFieldDriver;
 615     aMedRdFieldDriver.setFileName(filename_rd);
 616     aField_2->read( aMedRdFieldDriver );
 617   }
 618   //Test read(driverTypes driverType, const std::string & fileName);
 619   FIELD<double> * aField_3 = new FIELD<double>();
 620   aField_3->setName(fieldname_celldouble_rd);
 621   aField_3->read( MED_DRIVER, filename_rd);
 622
 623   ///////////////////
 624   //Test Write Part//
 625   ///////////////////
 626   //int IdDriver;
 627   MESH *aMesh = new MESH(MED_DRIVER,filename_rd,meshname);
 628   const SUPPORT *aSupport = aMesh->getSupportOnAll(MED_CELL);
 629   FIELD<int> *aFieldSupport;
 630   CPPUNIT_ASSERT_THROW(aFieldSupport =
 631                        new FIELD<int>(aSupport, MED_DRIVER,filename_rd,
 632                                       fieldname_nodeint_rd), MEDMEM::MEDEXCEPTION);
 633   aSupport = aMesh->getSupportOnAll(MED_NODE);
 634   CPPUNIT_ASSERT_NO_THROW(aFieldSupport =
 635                           new FIELD<int>(aSupport, MED_DRIVER, filename_rd,
 636                                          fieldname_nodeint_rd));
 637   aFieldSupport->removeReference();
 638
 639   //Test write(int index) method
 640   // Add drivers to FIELDs
 641   int IdDriver1 = aField_3->addDriver(MED_DRIVER,filename_wr,fieldname_celldouble_wr);
 642
 643   //Trying call write(int index) method with incorrect index
 644   CPPUNIT_ASSERT_THROW(aField_3->write(IdDriver1+1),MEDEXCEPTION);
 645
 646   //Write field to file
 647   aField_3->write(IdDriver1);
 648
 649   CPPUNIT_ASSERT_NO_THROW(aField_3->rmDriver(IdDriver1));
 650
 651   //Test write(const GENDRIVER &, MED_EN::med_mode_acces medMode=MED_EN::RDWR);
 652   {
 653     MED_FIELD_WRONLY_DRIVER<int> aMedWrFieldDriver;
 654     aMedWrFieldDriver.setFileName(filename_wr);
 655
 656     aField_3->setName(fieldname_nodeint_wr);
 657     aField_3->write(aMedWrFieldDriver);
 658     FIELD<double> aField_3_RD;
 659     aField_3_RD.setName(fieldname_nodeint_wr);
 660     aField_3_RD.read(MED_DRIVER,filename_wr);
 661   }
 662
 663   // Test write(driverTypes driverType, const std::string& filename
 664   //            MED_EN::med_mode_acces medMode=MED_EN::RDWR)
 665   aField_3->setName(fieldname_nodeint_wr1);
 666   // wrong mode
 667   CPPUNIT_ASSERT_THROW(aField_3->write(MED_DRIVER,filename_wr,MED_EN::RDONLY), MEDEXCEPTION);
 668   aField_3->write(MED_DRIVER,filename_wr);
 669   {
 670     FIELD<double> aField_3_RD;
 671     aField_3_RD.setName(fieldname_nodeint_wr1);
 672     aField_3_RD.read(MED_DRIVER,filename_wr);
 673   }
 674   aField_3->setName("fieldname_nodeint_wr1");
 675   aField_3->write(MED_DRIVER,filename_wr, MED_EN::WRONLY);// overwrite filename_wr
 676   {
 677     FIELD<double> aField_3_RD;
 678     aField_3_RD.setName(fieldname_nodeint_wr1);
 679     CPPUNIT_ASSERT_THROW(aField_3_RD.read(MED_DRIVER,filename_wr),MEDEXCEPTION);
 680   }
 681   //Test writeAppend(int index) method
 682   //Create a vtk file
 683   MESH * aMesh_1 = new MESH(MED_DRIVER,filename22_rd, meshname);
 684   aMesh_1->write(VTK_DRIVER, filenamevtk_wr);
 685   //Create a field
 686   FIELD<int> * aField_4 =
 687     new FIELD<int>(MED_DRIVER,filename22_rd,fieldname_nodeint_rd,-1,-1,aMesh_1);
 688   //Add Driver to a field
 689   int IdDriver2 = aField_4->addDriver(VTK_DRIVER, filenamevtk_wr ,fieldname_nodeint_wr);
 690   //Trying call writeAppend() method with incorrect index
 691   CPPUNIT_ASSERT_THROW(aField_4->writeAppend(IdDriver2+1,fieldname_nodeint_wr),MEDEXCEPTION);
 692
 693   CPPUNIT_ASSERT_NO_THROW(aField_4->writeAppend(IdDriver2, fieldname_nodeint_wr));
 694
 695   //Test writeAppend(const GENDRIVER &) method
 696   aField_4->setName(fieldname_nodeint_wr1);
 697   {
 698     VTK_FIELD_DRIVER<int> aVtkFieldDriver(filenamevtk_wr, aField_4);
 699     CPPUNIT_ASSERT_NO_THROW(aField_4->writeAppend(aVtkFieldDriver));
 700   }
 701
 702   //Delete objects
 703   aField_1->removeReference();
 704   aField_2->removeReference();
 705   aField_3->removeReference();
 706   aField_4->removeReference();
 707   aMesh->removeReference();
 708   aMesh_1->removeReference();
 709 }
 710
 711 void MEDMEMTest::testField()
 712 {
 713   SUPPORT *anEmptySupport=new SUPPORT;
 714   ////////////////////
 715   // TEST 1: FIELD_ //
 716   ////////////////////
 717   FIELD_ *aField_=new FIELD_;
 718
 719   // check set/get methods
 720   MED_EN::med_type_champ aValueType = MED_EN::MED_UNDEFINED_TYPE;
 721   MED_EN::medModeSwitch  anInterlace = MED_EN::MED_UNDEFINED_INTERLACE;
 722   checkField_(aField_, anEmptySupport, aValueType, anInterlace);
 723
 724   // copy constructor
 725   // This fails (Segmentation fault) if not set:
 726   // _componentsNames or _componentsDescriptions, or _MEDComponentsUnits
 727   FIELD_ *aField_copy1=new FIELD_(*aField_);
 728   compareField_(aField_, aField_copy1, /*isFIELD = */false, /*isValue = */false);
 729   aField_copy1->removeReference();
 730   // operator=
 731   FIELD_ *aField_copy2=new FIELD_;
 732   *aField_copy2 = *aField_;
 733   compareField_(aField_, aField_copy2,/*isFIELD = */false, /*isValue = */false);
 734   aField_copy2->removeReference();
 735   aField_->removeReference();
 736
 737   ////////////////////////
 738   // TEST 2: FIELD<int> //
 739   ////////////////////////
 740   FIELD<int> *aFieldInt=new FIELD<int>();
 741   checkField(aFieldInt, anEmptySupport);
 742   aFieldInt->removeReference();
 743   anEmptySupport->removeReference();
 744   ////////////////////////////////////////
 745   // TEST 3: FIELD<double, NoInterlace> //
 746   ////////////////////////////////////////
 747   MESH * aMesh = MEDMEMTest_createTestMesh();
 748   const GROUP* aGroup = aMesh->getGroup(MED_EN::MED_FACE, 1);
 749
 750   FIELD<double, NoInterlace> *aFieldDouble=new FIELD<double, NoInterlace>();
 751   checkField(aFieldDouble, aGroup);
 752   aFieldDouble->removeReference();
 753   //////////////////////////////////////////
 754   // TEST 4: FIELD<double, FullInterlace> //
 755   //////////////////////////////////////////
 756   FIELD<double> * aFieldOnGroup1 = createFieldOnGroup<double>(aMesh, aGroup, "Linear", "N");
 757   FIELD<double> * aFieldOnGroup2 = createFieldOnGroup<double>(aMesh, aGroup, "Quadratic", "N**2");
 758
 759   int nbVals = aFieldOnGroup1->getNumberOfValues();
 760   CPPUNIT_ASSERT(nbVals);
 761
 762   // numbers of elements in group,
 763   // they are needed in method FIELD::setValueIJ()
 764   const int *anElems = aGroup->getnumber()->getValue();
 765   double eucl1 = 0., eucl2 = 0.;
 766
 767   for (int i = 1; i <= nbVals; i++)
 768     {
 769       aFieldOnGroup1->setValueIJ(anElems[i-1], 1, (double)i);
 770       aFieldOnGroup1->setValueIJ(anElems[i-1], 2, (double)(-i));
 771
 772       aFieldOnGroup2->setValueIJ(anElems[i-1], 1, (double)i*i);
 773       aFieldOnGroup2->setValueIJ(anElems[i-1], 2, (double)(-i*i));
 774
 775       eucl1 += 2. * i * i;
 776       eucl2 += 2. * i * i * i * i;
 777     }
 778
 779   // out of bound (inexisting 33-th component of last element)
 780   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->setValueIJ(anElems[nbVals-1], 33, 10.), MEDEXCEPTION);
 781
 782   // normMax
 783   CPPUNIT_ASSERT_DOUBLES_EQUAL(nbVals, aFieldOnGroup1->normMax(), 0.000001);
 784   CPPUNIT_ASSERT_DOUBLES_EQUAL(nbVals*nbVals, aFieldOnGroup2->normMax(), 0.000001);
 785
 786   // norm2
 787   CPPUNIT_ASSERT_DOUBLES_EQUAL(std::sqrt(eucl1), aFieldOnGroup1->norm2(), 0.000001); // 10.4881
 788   CPPUNIT_ASSERT_DOUBLES_EQUAL(std::sqrt(eucl2), aFieldOnGroup2->norm2(), 0.000001); // 44.2493
 789
 790   // check getXXX methods
 791   CPPUNIT_ASSERT(!aFieldOnGroup1->getGaussPresence());
 792   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getArrayGauss(), MEDEXCEPTION);
 793   MEDMEM_ArrayInterface<double, FullInterlace,NoGauss>::Array * anArrayNoGauss =
 794     aFieldOnGroup1->getArrayNoGauss();
 795
 796   MEDMEM_Array_ * aMEDMEM_Array_ = aFieldOnGroup1->getArray();
 797   MEDMEM_ArrayInterface<double, FullInterlace,NoGauss>::Array * aMEDMEM_Array_conv =
 798     static_cast<MEDMEM_ArrayInterface<double, FullInterlace,NoGauss>::Array *>(aMEDMEM_Array_);
 799
 800   const double * aValues = aFieldOnGroup1->getValue();
 801
 802   // out of range
 803   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getColumn(3), MEDEXCEPTION);
 804   // cannot get column in FullInterlace
 805   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getColumn(1), MEDEXCEPTION);
 806
 807   for (int i = 1; i <= nbVals; i++)
 808     {
 809       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , aFieldOnGroup1->getValueIJK(anElems[i-1], 1, 1), 0.000001);
 810       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), aFieldOnGroup1->getValueIJK(anElems[i-1], 2, 1), 0.000001);
 811
 812       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , aValues[(i-1)*2 + 0], 0.000001);
 813       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), aValues[(i-1)*2 + 1], 0.000001);
 814
 815       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , anArrayNoGauss->getIJ(i, 1), 0.000001);
 816       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), anArrayNoGauss->getIJ(i, 2), 0.000001);
 817
 818       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , aMEDMEM_Array_conv->getIJ(i, 1), 0.000001);
 819       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), aMEDMEM_Array_conv->getIJ(i, 2), 0.000001);
 820
 821       const double* row_i = aFieldOnGroup1->getRow(anElems[i-1]);
 822       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , row_i[0], 0.000001);
 823       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), row_i[1], 0.000001);
 824
 825       double vals_i [2];
 826       aFieldOnGroup1->getValueOnElement(anElems[i-1], vals_i);
 827       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)i   , vals_i[0], 0.000001);
 828       CPPUNIT_ASSERT_DOUBLES_EQUAL((double)(-i), vals_i[1], 0.000001);
 829     }
 830   {
 831     // Test getValueOnPoints()
 832     // -----------------------
 833     const string file = getResourceFile("cube_hexa8.med");
 834     MESH mesh(MED_DRIVER, file, "CUBE_EN_HEXA8");
 835     double value[6];
 836     {
 837       FIELD<double> fieldcelldouble(MED_DRIVER, file, "fieldcelldouble", -1,-1, &mesh);
 838       fieldcelldouble.getSupport()->setMesh( &mesh );
 839       double point[6] =
 840         {
 841           0.25, 0.75, 0.25, // the 3-d cell
 842           1.0,  1.0,  1.0
 843         };// the 8-th cell
 844       fieldcelldouble.getValueOnPoint(point, value);
 845       const double tol = 1e-20;
 846       CPPUNIT_ASSERT_DOUBLES_EQUAL( 2.0, value[0], tol );
 847       CPPUNIT_ASSERT_DOUBLES_EQUAL( 3.0, value[1], tol );
 848       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, value[2], tol );
 849       fieldcelldouble.getValueOnPoints(2, point, value);
 850       CPPUNIT_ASSERT_DOUBLES_EQUAL( 2.0, value[0], tol );
 851       CPPUNIT_ASSERT_DOUBLES_EQUAL( 3.0, value[1], tol );
 852       CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, value[2], tol );
 853       CPPUNIT_ASSERT_DOUBLES_EQUAL( 2.0, value[3], tol );
 854       CPPUNIT_ASSERT_DOUBLES_EQUAL( 2.0, value[4], tol );
 855       CPPUNIT_ASSERT_DOUBLES_EQUAL( 3.0, value[5], tol );
 856     }
 857     {
 858       FIELD<double> fieldnodedouble(MED_DRIVER, file, "fieldnodedouble", -1,-1, &mesh); // t=0.0
 859       fieldnodedouble.getSupport()->setMesh( &mesh );
 860       double point[6] =
 861         {
 862           1.0, 1.0, 0.0, // the 9-th node
 863           1.0, 0.0, 1.0
 864         };// the 21-st node
 865       fieldnodedouble.getValueOnPoint(point, value);
 866       const double tol = 1e-20;
 867       CPPUNIT_ASSERT_DOUBLES_EQUAL( 5.0, value[0], tol );
 868       fieldnodedouble.getValueOnPoints(2, point, value);
 869       CPPUNIT_ASSERT_DOUBLES_EQUAL( 5.0, value[0], tol );
 870       CPPUNIT_ASSERT_DOUBLES_EQUAL( 4.0, value[1], tol );
 871
 872       point[0] = 1.1; // point out of mesh
 873       CPPUNIT_ASSERT_THROW(fieldnodedouble.getValueOnPoint(point, value), MEDEXCEPTION);
 874     }
 875   }
 876
 877   // modify all values of aFieldOnGroup2 by formula a*x + b (a = 2, b = 3)
 878   aFieldOnGroup2->applyLin(2., 3.);
 879   for (int i = 1; i <= nbVals; i++)
 880     {
 881       CPPUNIT_ASSERT_DOUBLES_EQUAL(3. + 2.*i*i, aFieldOnGroup2->getValueIJ(anElems[i-1], 1), 0.000001);
 882       CPPUNIT_ASSERT_DOUBLES_EQUAL(3. - 2.*i*i, aFieldOnGroup2->getValueIJ(anElems[i-1], 2), 0.000001);
 883     }
 884
 885   // apply function plus13() to aFieldOnGroup1
 886   aFieldOnGroup1->applyFunc<plus13>();
 887   for (int i = 1; i <= nbVals; i++)
 888     {
 889       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
 890       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
 891     }
 892
 893   // scalarProduct
 894   FIELD<double, FullInterlace> * aScalarProduct =
 895     FIELD<double, FullInterlace>::scalarProduct(*aFieldOnGroup1, *aFieldOnGroup2, /*deepCheck = */true);
 896   CPPUNIT_ASSERT_EQUAL(nbVals, aScalarProduct->getNumberOfValues());
 897   CPPUNIT_ASSERT_EQUAL(1, aScalarProduct->getNumberOfComponents());
 898   for (int i = 1; i <= nbVals; i++)
 899     {
 900       CPPUNIT_ASSERT_DOUBLES_EQUAL(78. + 4.*i*i*i, //(3. + 2.*i*i)*(13 + i) + (3. - 2.*i*i)*(13 - i)
 901                                    aScalarProduct->getValueIJ(anElems[i-1], 1), 0.000001);
 902     }
 903
 904   // fillFromAnalytic
 905   aFieldOnGroup2->fillFromAnalytic(proj2d);
 906
 907   double bary [3];
 908   double outp [2];
 909   const SUPPORT * aSupp = aFieldOnGroup2->getSupport();
 910   FIELD<double,FullInterlace> * barycenter = aMesh->getBarycenter(aSupp);
 911   for (int i = 1; i <= nbVals; i++)
 912     {
 913       bary[0] = barycenter->getValueIJ(anElems[i-1], 1);
 914       bary[1] = barycenter->getValueIJ(anElems[i-1], 2);
 915       bary[2] = barycenter->getValueIJ(anElems[i-1], 3);
 916
 917       proj2d(bary, outp);
 918
 919       //cout << "barycenter (" << bary[0] << ", " << bary[1] << ", " << bary[2] << ")" << endl;
 920       //cout << "proj2d     (" << outp[0] << ", " << outp[1] << ")" << endl;
 921
 922       //bary (-0.666667,  0.666667, 0.666667) -> outp ( 0.902369, 0.902369)
 923       //bary ( 0.666667, -0.666667, 0.666667) -> outp (-0.902369, 0.430964)
 924       //bary ( 0.      ,  0.      , 2.      ) -> outp ( 0.      , 2.      )
 925       //bary ( 0.      ,  0.      , 3.      ) -> outp ( 0.      , 3.      )
 926       //bary (-1.      ,  0.      , 2.5     ) -> outp ( 0.353553, 2.85355 )
 927
 928       //#ifdef ENABLE_FORCED_FAILURES
 929       // (BUG) in FIELD::fillFromAnalytic() in case of support, different from nodes:
 930       //       barycenterField in FullInterlace, but values extracted like from NoInterlace
 931       // TODO: fix MEDMEM_Field:3483 - code should depend on interlace
 932       CPPUNIT_ASSERT_DOUBLES_EQUAL(outp[0], aFieldOnGroup2->getValueIJ(anElems[i-1], 1), 0.000001);
 933       CPPUNIT_ASSERT_DOUBLES_EQUAL(outp[1], aFieldOnGroup2->getValueIJ(anElems[i-1], 2), 0.000001);
 934       //#endif
 935
 936       // currently it gives values, that are wrong:
 937       //aFieldOnGroup2 row1 ( 0.902369,  0.235702)
 938       //aFieldOnGroup2 row2 (-0.235702,  2.7643  )
 939       //aFieldOnGroup2 row3 (-0.235702, -1.2357  )
 940       //aFieldOnGroup2 row4 ( 1.7643  , -0.235702)
 941       //aFieldOnGroup2 row5 ( 0.235702,  2.7357  )
 942     }
 943
 944   // info about support (Group1)
 945   CPPUNIT_ASSERT(!aFieldOnGroup1->isOnAllElements()); // because we build Group1 so
 946   int nbTypes = aFieldOnGroup1->getNumberOfGeometricTypes();
 947   //CPPUNIT_ASSERT(nbTypes);
 948   CPPUNIT_ASSERT_EQUAL(2, nbTypes);
 949   const int * nbElemsInEachType = aFieldOnGroup1->getNumberOfElements();
 950   const MED_EN::medGeometryElement * aGeomTypes = aFieldOnGroup1->getGeometricTypes();
 951
 952   CPPUNIT_ASSERT_EQUAL(MED_EN::MED_TRIA3, aGeomTypes[0]);
 953   CPPUNIT_ASSERT_EQUAL(MED_EN::MED_QUAD4, aGeomTypes[1]);
 954
 955   // GAUSS
 956
 957   // now we have no gauss localization in aFieldOnGroup1
 958   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA3));
 959   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_QUAD4));
 960   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA6), MEDEXCEPTION);
 961   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(), MEDEXCEPTION);
 962
 963   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalization(MED_EN::MED_TRIA3), MEDEXCEPTION);
 964   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalizationPtr(MED_EN::MED_TRIA3), MEDEXCEPTION);
 965
 966   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNbGaussI(anElems[0]));
 967
 968   // set a gauss localization into aFieldOnGroup1
 969   double cooRef[6] =
 970     {
 971       1.,1., 2.,4., 3.,9.
 972     }; // xy xy xy
 973   double cooGauss[10] =
 974     {
 975       7.,7., 6.,6., 5.,5., 4.,3., 2.,1.
 976     }; // x1,y1  x2,y2  x3,y3  x4,y4  x5,y5
 977   double wg[5] =
 978     {
 979       1., 2., 3., 4., 5.
 980     };
 981   GAUSS_LOCALIZATION<> gl1 ("GL1", MED_EN::MED_TRIA3, /*nGauss*/5, cooRef, cooGauss, wg);
 982
 983   aFieldOnGroup1->setGaussLocalization(MED_EN::MED_TRIA3, gl1);
 984
 985   // now we have a gauss localization for MED_TRIA3 type
 986   CPPUNIT_ASSERT_EQUAL(5, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA3));
 987   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_QUAD4));
 988   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(MED_EN::MED_TRIA6), MEDEXCEPTION);
 989   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getNumberOfGaussPoints(), MEDEXCEPTION);
 990
 991   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalization(MED_EN::MED_QUAD4), MEDEXCEPTION);
 992   CPPUNIT_ASSERT_THROW(aFieldOnGroup1->getGaussLocalizationPtr(MED_EN::MED_QUAD4), MEDEXCEPTION);
 993
 994   GAUSS_LOCALIZATION<> gl1Back = aFieldOnGroup1->getGaussLocalization(MED_EN::MED_TRIA3);
 995   const GAUSS_LOCALIZATION<> * gl1BackPtr = aFieldOnGroup1->getGaussLocalizationPtr(MED_EN::MED_TRIA3);
 996
 997   CPPUNIT_ASSERT(gl1 == gl1Back);
 998   CPPUNIT_ASSERT(gl1 == *gl1BackPtr);
 999
1000   CPPUNIT_ASSERT_EQUAL(1, aFieldOnGroup1->getNbGaussI(anElems[0]));
1001
1002   // sub-support of Group1 on one (first) geometric type
1003   SUPPORT * aSubSupport1 = new SUPPORT();
1004   aSubSupport1->setMesh( aMesh );
1005   aSubSupport1->setName( "Sub-Support 1 of Group1" );
1006   aSubSupport1->setEntity( MED_EN::MED_FACE );
1007
1008   int nbTypes1 = 1;
1009   int nbElemsInEachType1[1];
1010   nbElemsInEachType1[0] = nbElemsInEachType[0];
1011   int nbElems1 = nbElemsInEachType1[0];
1012   MED_EN::medGeometryElement aGeomTypes1[1];
1013   aGeomTypes1[0] = aGeomTypes[0];
1014   int * anElems1 = new int[nbElems1];
1015   for (int i = 0; i < nbElems1; i++)
1016     {
1017       anElems1[i] = anElems[i];
1018     }
1019
1020   aSubSupport1->setpartial("Support for sub-field 1 on one type of elements",
1021                            nbTypes1, nbElems1, aGeomTypes1, nbElemsInEachType1, anElems1);
1022
1023   // extract sub-field on aSubSupport1
1024   FIELD<double, FullInterlace> * aSubField1 = aFieldOnGroup1->extract(aSubSupport1);
1025   CPPUNIT_ASSERT_EQUAL(nbElems1 * /*NumberOfComponents = */2, aSubField1->getValueLength());
1026
1027   // aSubField1:
1028   // elt\comp |  1 |  2
1029   //--------------------
1030   //  1       | 14 | 12
1031   //  2       | 15 | 11
1032
1033   // check normL2() and normL1()
1034   FIELD<double>* anAreaField = aMesh->getArea(aSubSupport1);
1035   double area1 = anAreaField->getValueIJ(anElems1[0], 1);
1036   double area2 = anAreaField->getValueIJ(anElems1[1], 1);
1037   CPPUNIT_ASSERT_DOUBLES_EQUAL(2.44949, area1, 0.00001);
1038   CPPUNIT_ASSERT_DOUBLES_EQUAL(2.44949, area2, 0.00001);
1039
1040   CPPUNIT_ASSERT_DOUBLES_EQUAL(210.5, aSubField1->normL2(1), 0.00001); // (14*14 + 15*15)/2
1041   CPPUNIT_ASSERT_DOUBLES_EQUAL(132.5, aSubField1->normL2(2), 0.00001); // (12*12 + 11*11)/2
1042   CPPUNIT_ASSERT_DOUBLES_EQUAL(343.0, aSubField1->normL2() , 0.00001); // 210.5 + 132.5
1043
1044   CPPUNIT_ASSERT_DOUBLES_EQUAL(14.5, aSubField1->normL1(1), 0.00001); // (14 + 15)/2
1045   CPPUNIT_ASSERT_DOUBLES_EQUAL(11.5, aSubField1->normL1(2), 0.00001); // (12 + 11)/2
1046   CPPUNIT_ASSERT_DOUBLES_EQUAL(26.0, aSubField1->normL1() , 0.00001); // 14.5 + 11.5
1047
1048   double aNewArea [2] =
1049     {
1050       1., 0.
1051     }; // only first element will be taken into account
1052   anAreaField->setValue(aNewArea);
1053
1054   CPPUNIT_ASSERT_DOUBLES_EQUAL(196.0, aSubField1->normL2(1, anAreaField), 0.00001); // 14*14/1
1055   CPPUNIT_ASSERT_DOUBLES_EQUAL(144.0, aSubField1->normL2(2, anAreaField), 0.00001); // 12*12/1
1056   CPPUNIT_ASSERT_DOUBLES_EQUAL(340.0, aSubField1->normL2(anAreaField) , 0.00001); // 196 + 144
1057
1058   CPPUNIT_ASSERT_DOUBLES_EQUAL(14.0, aSubField1->normL1(1, anAreaField), 0.00001); // 14/1
1059   CPPUNIT_ASSERT_DOUBLES_EQUAL(12.0, aSubField1->normL1(2, anAreaField), 0.00001); // 12/1
1060   CPPUNIT_ASSERT_DOUBLES_EQUAL(26.0, aSubField1->normL1(anAreaField) , 0.00001); // 14 + 12
1061
1062   // check normL2() on nodal field (issue 0020120)
1063   {
1064     // read nodal field from pointe.med
1065     string filename  = getResourceFile("pointe.med");
1066     string fieldname = "fieldnodedouble";
1067     string meshname  = "maa1";
1068     FIELD<double> *nodalField=new FIELD<double>( MED_DRIVER, filename, fieldname);
1069     MESH *mesh=new MESH( MED_DRIVER, filename, meshname);
1070     nodalField->getSupport()->setMesh( mesh );
1071
1072     // make a field on the nodes of first cell only
1073     SUPPORT *oneCellNodesSup=new SUPPORT;
1074     oneCellNodesSup->setMesh(mesh);
1075     oneCellNodesSup->setName("Sub-Support of nodes of 1 cell");
1076     oneCellNodesSup->setEntity(MED_NODE);
1077     int NumberOfElements[] =
1078       {
1079         mesh->getTypes(MED_CELL)[0]%100
1080       };
1081     medGeometryElement GeometricType[] =
1082       {
1083         MED_POINT1
1084       };
1085     oneCellNodesSup->setpartial("Support for sub-field of one cell nodes",
1086                                 /*NumberOfGeometricType=*/1,
1087                                 /*TotalNumberOfElements=*/ *NumberOfElements,
1088                                 GeometricType,
1089                                 NumberOfElements,
1090                                 /*NumberValue=*/ mesh->getConnectivity(MED_NODAL,MED_CELL,MED_ALL_ELEMENTS ));
1091     FIELD<double, FullInterlace> * oneCellNodesField = nodalField->extract( oneCellNodesSup );
1092     oneCellNodesSup->removeReference();
1093     // compute normL2 by avarage nodal value on the cell
1094
1095     const SUPPORT *allCellsSupport=mesh->getSupportOnAll( MED_CELL );
1096     FIELD<double>* volumeField = mesh->getVolume(allCellsSupport);
1097     // mdump output:
1098     // - Mailles de type MED_TETRA4 :
1099     //  [     1 ] :          1          2          3          6
1100     // - Valeurs :
1101     // | 1.000000 | 1.000000 | 1.000000 | 2.000000 | 2.000000 | 2.000000 | ...
1102     const double cellVal = ( 1.000000 + 1.000000 + 1.000000 + 2.000000 ) / 4;
1103     const double vol = abs( volumeField->getValueIJ( 1, 1 ));
1104     const double norm = cellVal*cellVal*vol/vol; // v*v*vol/totVol;
1105     CPPUNIT_ASSERT_DOUBLES_EQUAL(norm, oneCellNodesField->normL2(), 0.000001);
1106     CPPUNIT_ASSERT_DOUBLES_EQUAL(norm, oneCellNodesField->normL2(1), 0.000001);
1107     CPPUNIT_ASSERT_DOUBLES_EQUAL(norm, oneCellNodesField->normL2(volumeField), 0.000001);
1108     CPPUNIT_ASSERT_DOUBLES_EQUAL(norm, oneCellNodesField->normL2(1, volumeField), 0.000001);
1109
1110     CPPUNIT_ASSERT_THROW(oneCellNodesField->normL1(), MEDEXCEPTION);
1111     CPPUNIT_ASSERT_THROW(oneCellNodesField->normL1(1), MEDEXCEPTION);
1112     CPPUNIT_ASSERT_THROW(oneCellNodesField->normL2(nodalField), MEDEXCEPTION);
1113     CPPUNIT_ASSERT_THROW(oneCellNodesField->normL2(anAreaField), MEDEXCEPTION);
1114     CPPUNIT_ASSERT_THROW(oneCellNodesField->normL2(aSubField1), MEDEXCEPTION);
1115     CPPUNIT_ASSERT_THROW(oneCellNodesField->normL2(aFieldOnGroup1), MEDEXCEPTION);
1116     nodalField->removeReference();
1117     volumeField->removeReference();
1118     oneCellNodesField->removeReference();
1119     mesh->removeReference();
1120   }
1121
1122   // double integral(SUPPORT*) - mantis issue 0020460:
1123   // [CEA 352] Integral calculus on all field or on a field subarea (groupe or family)
1124   {
1125     // make 2D grid 2x2 with steps along axis 1.0, 2.0
1126     const int dim = 2;
1127     double coord[3] =
1128       {
1129         0.0, 1.0, 3.0
1130       };
1131     vector< vector<double> > xyz_array( dim, vector<double>( coord, coord + 3 ));
1132     vector<string> coord_name( dim, "coord_name");
1133     vector<string> coord_unit( dim, "m");
1134     MESH *mesh= const_cast<MESH*>( GRID( xyz_array, coord_name, coord_unit ).convertInMESH());
1135
1136     // make supports on the grid
1137     const SUPPORT *supOnAll=mesh->getSupportOnAll(MED_CELL);
1138     SUPPORT *sup123=new SUPPORT;
1139     sup123->setMesh(mesh);
1140     sup123->setName("sup123");
1141     int nbEl123[] =
1142       {
1143         3
1144       };
1145     int elems123[] =
1146       {
1147         1,2,3
1148       };
1149     SUPPORT *sup12=new SUPPORT;
1150     sup12->setMesh(mesh);
1151     sup12->setName("sup12");
1152     int nbEl12 [] =
1153       {
1154         2
1155       };
1156     int elems12 [] =
1157       {
1158         1,2
1159       };
1160     SUPPORT *sup34=new SUPPORT;
1161     sup34->setMesh(mesh);
1162     sup34->setName("sup34");
1163     int nbEl34 [] =
1164       {
1165         2
1166       };
1167     int elems34 [] =
1168       {
1169         3,4
1170       };
1171     const int nbGeomTypes = 1;
1172     const medGeometryElement * geomType = mesh->getTypes(MED_EN::MED_CELL);
1173     mesh->removeReference();
1174     sup123->setpartial("test", nbGeomTypes, *nbEl123, geomType, nbEl123, elems123 );
1175     sup12->setpartial("test", nbGeomTypes, *nbEl12 , geomType, nbEl12 , elems12 );
1176     sup34->setpartial("test", nbGeomTypes, *nbEl34 , geomType, nbEl34 , elems34 );
1177
1178     // make vectorial fields with values of i-th elem { i, i*10, i*100 }
1179     const int nbComp = 3, nbElems = 4;
1180     const int mult[nbComp] =
1181       {
1182         1, 10, 100
1183       };
1184     FIELD<int,NoInterlaceByType> *fAllNoTy=new FIELD<int,NoInterlaceByType>(supOnAll, nbComp), *f12NoTy=new FIELD<int,NoInterlaceByType>(sup12, nbComp);
1185     FIELD<int,NoInterlace>       *fAllNo=new FIELD<int,NoInterlace>(supOnAll, nbComp), *f12No=new FIELD<int,NoInterlace>(sup12, nbComp);
1186     FIELD<int,FullInterlace>     *fAllFull=new FIELD<int,FullInterlace>(supOnAll, nbComp), *f12Full=new FIELD<int,FullInterlace>(sup12, nbComp);
1187     int i, j;
1188     for ( i = 1; i <= nbElems; ++i )
1189       for ( j = 1; j <= nbComp; ++j )
1190         {
1191           fAllFull->setValueIJ( i, j, i * mult[j-1]);
1192           fAllNoTy->setValueIJ( i, j, i * mult[j-1]);
1193           fAllNo  ->setValueIJ( i, j, i * mult[j-1]);
1194           if ( i < 3 )
1195             {
1196               f12Full->setValueIJ( i, j, i * mult[j-1]);
1197               f12NoTy->setValueIJ( i, j, i * mult[j-1]);
1198               f12No  ->setValueIJ( i, j, i * mult[j-1]);
1199             }
1200         }
1201     // Test
1202     double preci = 1e-18, integral;
1203     // Integral = SUM( area * (i*1 + i*10 + i*100)) == 111 * SUM( area * i )
1204     // elem area: { 1, 2, 2, 4 }
1205     integral = 111*( 1*1 + 2*2 + 2*3 + 4*4 );
1206     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllNoTy->integral(), preci );
1207     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllNo  ->integral(), preci );
1208     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllFull->integral(), preci );
1209     integral = 111*( 1*1 + 2*2 );
1210     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12NoTy->integral(), preci );
1211     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12No  ->integral(), preci );
1212     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12Full->integral(), preci );
1213
1214     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllNoTy->integral(sup12), preci );
1215     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllNo  ->integral(sup12), preci );
1216     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllFull->integral(sup12), preci );
1217
1218     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12NoTy->integral(sup12), preci );
1219     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12No  ->integral(sup12), preci );
1220     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12Full->integral(sup12), preci );
1221     sup12->removeReference();
1222     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12NoTy->integral(sup123), preci );
1223     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12No  ->integral(sup123), preci );
1224     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12Full->integral(sup123), preci );
1225     integral = 111*( 1*1 + 2*2 + 2*3 );
1226     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllNoTy->integral(sup123), preci );
1227     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllNo  ->integral(sup123), preci );
1228     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, fAllFull->integral(sup123), preci );
1229     fAllNoTy->removeReference();
1230     fAllNo->removeReference();
1231     sup123->removeReference();
1232     fAllFull->removeReference();
1233     integral = 0;
1234     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12NoTy->integral(sup34), preci );
1235     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12No  ->integral(sup34), preci );
1236     CPPUNIT_ASSERT_DOUBLES_EQUAL( integral, f12Full->integral(sup34), preci );
1237     sup34->removeReference();
1238     f12NoTy->removeReference();
1239     f12No->removeReference();
1240     f12Full->removeReference();
1241   }
1242
1243   // applyPow
1244   aSubField1->applyPow(2.);
1245   CPPUNIT_ASSERT_DOUBLES_EQUAL(196., aSubField1->getValueIJ(anElems1[0], 1), 0.000001); // 14*14
1246   CPPUNIT_ASSERT_DOUBLES_EQUAL(144., aSubField1->getValueIJ(anElems1[0], 2), 0.000001); // 12*12
1247   CPPUNIT_ASSERT_DOUBLES_EQUAL(225., aSubField1->getValueIJ(anElems1[1], 1), 0.000001); // 15*15
1248   CPPUNIT_ASSERT_DOUBLES_EQUAL(121., aSubField1->getValueIJ(anElems1[1], 2), 0.000001); // 11*11
1249
1250   // setArray (NoGauss)
1251   MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array * aNewArrayNoGauss =
1252     new MEDMEM_ArrayInterface<double,FullInterlace,NoGauss>::Array(/*dim*/2, /*nbelem*/2);
1253   aNewArrayNoGauss->setIJ(1, 1, 4.);
1254   aNewArrayNoGauss->setIJ(1, 2, 2.);
1255   aNewArrayNoGauss->setIJ(2, 1, 5.);
1256   aNewArrayNoGauss->setIJ(2, 2, 1.);
1257   aSubField1->setArray(aNewArrayNoGauss);
1258   // no need to delete aNewArrayNoGauss, because it will be deleted
1259   // in destructor or in deallocValue() method of aSubField1
1260
1261   CPPUNIT_ASSERT_DOUBLES_EQUAL(4., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
1262   CPPUNIT_ASSERT_DOUBLES_EQUAL(2., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
1263   CPPUNIT_ASSERT_DOUBLES_EQUAL(5., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
1264   CPPUNIT_ASSERT_DOUBLES_EQUAL(1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
1265
1266   // setRow
1267   double row[2] =
1268     {
1269       -1., -3.
1270     };
1271   aSubField1->setRow(anElems1[0], row);
1272   CPPUNIT_ASSERT_DOUBLES_EQUAL(-1., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
1273   CPPUNIT_ASSERT_DOUBLES_EQUAL(-3., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
1274   CPPUNIT_ASSERT_DOUBLES_EQUAL( 5., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
1275   CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
1276   // out of range
1277   CPPUNIT_ASSERT_THROW(aSubField1->setRow(3, row), MEDEXCEPTION);
1278
1279   // setColumn
1280   double col[2] =
1281     {
1282       -7., -9.
1283     };
1284   aSubField1->setColumn(1, col);
1285   CPPUNIT_ASSERT_DOUBLES_EQUAL(-7., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
1286   CPPUNIT_ASSERT_DOUBLES_EQUAL(-3., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
1287   CPPUNIT_ASSERT_DOUBLES_EQUAL(-9., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
1288   CPPUNIT_ASSERT_DOUBLES_EQUAL( 1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
1289   // out of range
1290   CPPUNIT_ASSERT_THROW(aSubField1->setColumn(3, col), MEDEXCEPTION);
1291
1292   // setArray (Gauss)
1293   {
1294     int nbelgeoc[2] =
1295       {
1296         1, 3
1297       }; // 3 - 1 = two elements for the first (and the only) type
1298     int nbgaussgeo[2] =
1299       {
1300         0, 1
1301       }; // one gauss point per each element
1302     MEDMEM_ArrayInterface<double,FullInterlace,Gauss>::Array * aNewArrayGauss =
1303       new MEDMEM_ArrayInterface<double,FullInterlace,Gauss>::Array
1304       (/*dim*/2, /*nbelem*/2, /*nbtypegeo*/1, /*nbelgeoc*/nbelgeoc, /*nbgaussgeo*/nbgaussgeo);
1305
1306     aNewArrayGauss->setIJ(1, 1, -4.);
1307     aNewArrayGauss->setIJ(1, 2, -2.);
1308     aNewArrayGauss->setIJ(2, 1, -5.);
1309     aNewArrayGauss->setIJ(2, 2, -1.);
1310
1311     aNewArrayGauss->setIJK(1, 1, 1, -4.);
1312     aNewArrayGauss->setIJK(1, 2, 1, -2.);
1313     aNewArrayGauss->setIJK(2, 1, 1, -5.);
1314     aNewArrayGauss->setIJK(2, 2, 1, -1.);
1315
1316     aSubField1->setArray(aNewArrayGauss);
1317     // no need to delete aNewArrayGauss, because it will be deleted
1318     // in destructor or in deallocValue() method of aSubField1
1319
1320     CPPUNIT_ASSERT_DOUBLES_EQUAL(-4., aSubField1->getValueIJ(anElems1[0], 1), 0.000001);
1321     CPPUNIT_ASSERT_DOUBLES_EQUAL(-2., aSubField1->getValueIJ(anElems1[0], 2), 0.000001);
1322     CPPUNIT_ASSERT_DOUBLES_EQUAL(-5., aSubField1->getValueIJ(anElems1[1], 1), 0.000001);
1323     CPPUNIT_ASSERT_DOUBLES_EQUAL(-1., aSubField1->getValueIJ(anElems1[1], 2), 0.000001);
1324
1325     CPPUNIT_ASSERT_DOUBLES_EQUAL(-4., aSubField1->getValueIJK(anElems1[0], 1, 1), 0.000001);
1326     CPPUNIT_ASSERT_DOUBLES_EQUAL(-2., aSubField1->getValueIJK(anElems1[0], 2, 1), 0.000001);
1327     CPPUNIT_ASSERT_DOUBLES_EQUAL(-5., aSubField1->getValueIJK(anElems1[1], 1, 1), 0.000001);
1328     CPPUNIT_ASSERT_DOUBLES_EQUAL(-1., aSubField1->getValueIJK(anElems1[1], 2, 1), 0.000001);
1329   }
1330
1331   // alloc/dealloc; compatibility of new size with support
1332   try
1333     {
1334       aSubField1->deallocValue();
1335       aSubField1->allocValue(/*NumberOfComponents*/2, /*LengthValue*/5);
1336       //#ifdef ENABLE_FAULTS
1337       // (BUG) No compatibility between Support and allocated value
1338       //aSubField1->normL1();
1339       CPPUNIT_ASSERT_THROW(aSubField1->normL1(),MEDEXCEPTION);
1340       //#endif
1341       //#ifdef ENABLE_FORCED_FAILURES
1342       // TODO: FIX to throw an EXCEPTION
1343       //    CPPUNIT_FAIL("Error: no compatibility between Support and allocated value");
1344       //#endif
1345     }
1346   catch (MEDEXCEPTION & ex)
1347     {
1348       // normal behaviour
1349     }
1350   catch (...)
1351     {
1352       CPPUNIT_FAIL("Error: no compatibility between Support and allocated value");
1353     }
1354
1355   // check that aFieldOnGroup1 is not changed after aSubField1 modifications
1356   for (int i = 1; i <= nbVals; i++)
1357     {
1358       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
1359       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
1360     }
1361
1362   // reset aFieldOnGroup2 values for simple control of operators results
1363   for (int i = 1; i <= nbVals; i++)
1364     {
1365       aFieldOnGroup2->setValueIJ(anElems[i-1], 1, i*i);
1366       aFieldOnGroup2->setValueIJ(anElems[i-1], 2, -i*i);
1367     }
1368
1369   int len = aFieldOnGroup1->getValueLength();
1370   const double * val1 = aFieldOnGroup1->getValue();
1371   const double * val2 = aFieldOnGroup2->getValue();
1372   const double * val_res;
1373
1374   // operators and add, sub, mul, div
1375
1376   // +
1377   FIELD<double> *aSum = *aFieldOnGroup1 + *aFieldOnGroup2;
1378   aSum->setName(aFieldOnGroup1->getName());
1379   aSum->setDescription(aFieldOnGroup1->getDescription());
1380   compareField_(aFieldOnGroup1, aSum, true, true);
1381   val_res = aSum->getValue();
1382   for (int i = 0; i < len; i++)
1383     {
1384       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] + val2[i], val_res[i], 0.000001);
1385     }
1386   aSum->removeReference();
1387   // -
1388   FIELD<double> *aDifference = *aFieldOnGroup1 - *aFieldOnGroup2;
1389   aDifference->setName(aFieldOnGroup1->getName());
1390   aDifference->setDescription(aFieldOnGroup1->getDescription());
1391   compareField_(aFieldOnGroup1, aDifference, true, true);
1392   val_res = aDifference->getValue();
1393   for (int i = 0; i < len; i++)
1394     {
1395       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] - val2[i], val_res[i], 0.000001);
1396     }
1397   aDifference->removeReference();
1398   // - (unary)
1399   FIELD<double> *aNegative = - *aFieldOnGroup1;
1400   aNegative->setName(aFieldOnGroup1->getName());
1401   aNegative->setDescription(aFieldOnGroup1->getDescription());
1402   compareField_(aFieldOnGroup1, aNegative, true, true);
1403   val_res = aNegative->getValue();
1404   for (int i = 0; i < len; i++)
1405     {
1406       CPPUNIT_ASSERT_DOUBLES_EQUAL(- val1[i], val_res[i], 0.000001);
1407     }
1408   aNegative->removeReference();
1409   // *
1410   FIELD<double> *aProduct = (*aFieldOnGroup1) * (*aFieldOnGroup2);
1411   aProduct->setName(aFieldOnGroup1->getName());
1412   aProduct->setDescription(aFieldOnGroup1->getDescription());
1413   compareField_(aFieldOnGroup1, aProduct, true, true);
1414   val_res = aProduct->getValue();
1415   for (int i = 0; i < len; i++)
1416     {
1417       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] * val2[i], val_res[i], 0.000001);
1418     }
1419   aProduct->removeReference();
1420   // /
1421   FIELD<double> *aQuotient = *aFieldOnGroup1 / *aFieldOnGroup2;
1422   aQuotient->setName(aFieldOnGroup1->getName());
1423   aQuotient->setDescription(aFieldOnGroup1->getDescription());
1424   compareField_(aFieldOnGroup1, aQuotient, true, true);
1425   val_res = aQuotient->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   aQuotient->removeReference();
1431   double val22 = aFieldOnGroup2->getValueIJ(anElems[2], 2);
1432   aFieldOnGroup2->setValueIJ(anElems[2], 2, 0.);
1433
1434   CPPUNIT_ASSERT_THROW((*aFieldOnGroup1 / *aFieldOnGroup2), MEDEXCEPTION);
1435   //#ifdef ENABLE_FORCED_FAILURES
1436   // (BUG) is it up to user to control validity of data to avoid division on zero?
1437   // YES: USER SHOULD CARE OF IT
1438   //CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 /= *aFieldOnGroup2, MEDEXCEPTION);
1439   //#endif
1440   CPPUNIT_ASSERT_THROW(FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup2)->removeReference(), MEDEXCEPTION);
1441   CPPUNIT_ASSERT_THROW(FIELD<double>::divDeep(*aFieldOnGroup1, *aFieldOnGroup2)->removeReference(), MEDEXCEPTION);
1442
1443   // restore value
1444   aFieldOnGroup2->setValueIJ(anElems[2], 2, val22);
1445
1446   // restore values
1447   for (int i = 1; i <= nbVals; i++)
1448     {
1449       aFieldOnGroup1->setValueIJ(anElems[i-1], 1, 13 + i);
1450       aFieldOnGroup1->setValueIJ(anElems[i-1], 2, 13 - i);
1451     }
1452
1453   // static methods
1454   FIELD<double> * aPtr;
1455
1456   // add
1457   aPtr = FIELD<double>::add(*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   aPtr->removeReference();
1467
1468   // sub
1469   aPtr = FIELD<double>::sub(*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   aPtr->removeReference();
1479
1480   // mul
1481   aPtr = FIELD<double>::mul(*aFieldOnGroup1, *aFieldOnGroup2);
1482   aPtr->setName(aFieldOnGroup1->getName());
1483   aPtr->setDescription(aFieldOnGroup1->getDescription());
1484   compareField_(aFieldOnGroup1, aPtr, true, true);
1485   val_res = aPtr->getValue();
1486   for (int i = 0; i < len; i++)
1487     {
1488       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] * val2[i], val_res[i], 0.000001);
1489     }
1490   aPtr->removeReference();
1491
1492   // div
1493   aPtr = FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup2);
1494   aPtr->setName(aFieldOnGroup1->getName());
1495   aPtr->setDescription(aFieldOnGroup1->getDescription());
1496   compareField_(aFieldOnGroup1, aPtr, true, true);
1497   val_res = aPtr->getValue();
1498   for (int i = 0; i < len; i++)
1499     {
1500       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] / val2[i], val_res[i], 0.000001);
1501     }
1502   aPtr->removeReference();
1503
1504   // addDeep
1505   aPtr = FIELD<double>::addDeep(*aFieldOnGroup1, *aFieldOnGroup2);
1506   aPtr->setName(aFieldOnGroup1->getName());
1507   aPtr->setDescription(aFieldOnGroup1->getDescription());
1508   compareField_(aFieldOnGroup1, aPtr, true, true);
1509   val_res = aPtr->getValue();
1510   for (int i = 0; i < len; i++)
1511     {
1512       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] + val2[i], val_res[i], 0.000001);
1513     }
1514   aPtr->removeReference();
1515
1516   // subDeep
1517   aPtr = FIELD<double>::subDeep(*aFieldOnGroup1, *aFieldOnGroup2);
1518   aPtr->setName(aFieldOnGroup1->getName());
1519   aPtr->setDescription(aFieldOnGroup1->getDescription());
1520   compareField_(aFieldOnGroup1, aPtr, true, true);
1521   val_res = aPtr->getValue();
1522   for (int i = 0; i < len; i++)
1523     {
1524       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] - val2[i], val_res[i], 0.000001);
1525     }
1526   aPtr->removeReference();
1527
1528   // mulDeep
1529   aPtr = FIELD<double>::mulDeep(*aFieldOnGroup1, *aFieldOnGroup2);
1530   aPtr->setName(aFieldOnGroup1->getName());
1531   aPtr->setDescription(aFieldOnGroup1->getDescription());
1532   compareField_(aFieldOnGroup1, aPtr, true, true);
1533   val_res = aPtr->getValue();
1534   for (int i = 0; i < len; i++)
1535     {
1536       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] * val2[i], val_res[i], 0.000001);
1537     }
1538   aPtr->removeReference();
1539
1540   // divDeep
1541   aPtr = FIELD<double>::divDeep(*aFieldOnGroup1, *aFieldOnGroup2);
1542   aPtr->setName(aFieldOnGroup1->getName());
1543   aPtr->setDescription(aFieldOnGroup1->getDescription());
1544   compareField_(aFieldOnGroup1, aPtr, true, true);
1545   val_res = aPtr->getValue();
1546   for (int i = 0; i < len; i++)
1547     {
1548       CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] / val2[i], val_res[i], 0.000001);
1549     }
1550   aPtr->removeReference();
1551
1552   // +=
1553   *aFieldOnGroup1 += *aFieldOnGroup2;
1554   for (int i = 1; i <= nbVals; i++)
1555     {
1556       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i + i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
1557       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i - i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
1558     }
1559
1560   // -=
1561   *aFieldOnGroup1 -= *aFieldOnGroup2;
1562   for (int i = 1; i <= nbVals; i++)
1563     {
1564       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
1565       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
1566     }
1567
1568   // *=
1569   *aFieldOnGroup1 *= *aFieldOnGroup2;
1570   for (int i = 1; i <= nbVals; i++)
1571     {
1572       CPPUNIT_ASSERT_DOUBLES_EQUAL( (13 + i)*i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
1573       CPPUNIT_ASSERT_DOUBLES_EQUAL(-(13 - i)*i*i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
1574     }
1575
1576   // /=
1577   *aFieldOnGroup1 /= *aFieldOnGroup2;
1578   for (int i = 1; i <= nbVals; i++)
1579     {
1580       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 + i, aFieldOnGroup1->getValueIJ(anElems[i-1], 1), 0.000001);
1581       CPPUNIT_ASSERT_DOUBLES_EQUAL(13 - i, aFieldOnGroup1->getValueIJ(anElems[i-1], 2), 0.000001);
1582     }
1583
1584   // check case of different operands: support
1585   MESH * aMeshOneMore = MEDMEMTest_createTestMesh();
1586   const GROUP* aGroupOneMore = aMeshOneMore->getGroup(MED_EN::MED_FACE, 1);
1587   FIELD<double> * aFieldOnGroup3 =
1588     createFieldOnGroup<double>(aMeshOneMore, aGroupOneMore, "Test_Diff_Mesh", "test");
1589   for (int i = 1; i <= nbVals; i++)
1590     {
1591       aFieldOnGroup3->setValueIJ(anElems[i-1], 1, 2*i);
1592       aFieldOnGroup3->setValueIJ(anElems[i-1], 2, 3*i);
1593     }
1594   const double * val3 = aFieldOnGroup3->getValue();
1595
1596   //CPPUNIT_ASSERT_NO_THROW();
1597   try
1598     {
1599       aPtr = FIELD<double>::addDeep(*aFieldOnGroup1, *aFieldOnGroup3);
1600       aPtr->setName(aFieldOnGroup1->getName());
1601       aPtr->setDescription(aFieldOnGroup1->getDescription());
1602       compareField_(aFieldOnGroup1, aPtr, true, true);
1603       val_res = aPtr->getValue();
1604       for (int i = 0; i < len; i++)
1605         {
1606           CPPUNIT_ASSERT_DOUBLES_EQUAL(val1[i] + val3[i], val_res[i], 0.000001);
1607         }
1608       aPtr->removeReference();
1609
1610       aPtr = FIELD<double>::subDeep(*aFieldOnGroup1, *aFieldOnGroup3);
1611       aPtr->removeReference();
1612       aPtr = FIELD<double>::mulDeep(*aFieldOnGroup1, *aFieldOnGroup3);
1613       aPtr->removeReference();
1614       aPtr = FIELD<double>::divDeep(*aFieldOnGroup1, *aFieldOnGroup3);
1615       aPtr->removeReference();
1616     }
1617   catch (MEDEXCEPTION & ex)
1618     {
1619       CPPUNIT_FAIL(ex.what());
1620     }
1621   catch (...)
1622     {
1623       CPPUNIT_FAIL("Unknown exception in FIELD::xxxDeep()");
1624     }
1625
1626   CPPUNIT_ASSERT_THROW(FIELD<double>::add(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
1627   CPPUNIT_ASSERT_THROW(FIELD<double>::sub(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
1628   CPPUNIT_ASSERT_THROW(FIELD<double>::mul(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
1629   CPPUNIT_ASSERT_THROW(FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup3), MEDEXCEPTION);
1630
1631   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 + *aFieldOnGroup3, MEDEXCEPTION);
1632   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 - *aFieldOnGroup3, MEDEXCEPTION);
1633   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 * *aFieldOnGroup3, MEDEXCEPTION);
1634   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 / *aFieldOnGroup3, MEDEXCEPTION);
1635
1636   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 += *aFieldOnGroup3, MEDEXCEPTION);
1637   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 -= *aFieldOnGroup3, MEDEXCEPTION);
1638   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 *= *aFieldOnGroup3, MEDEXCEPTION);
1639   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 /= *aFieldOnGroup3, MEDEXCEPTION);
1640
1641   // check case of different operands: MEDComponentsUnits
1642   aFieldOnGroup1->setMEDComponentUnit(1, "unit3");
1643
1644   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 + *aFieldOnGroup2, MEDEXCEPTION);
1645   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 - *aFieldOnGroup2, MEDEXCEPTION);
1646   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 += *aFieldOnGroup2, MEDEXCEPTION);
1647   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 -= *aFieldOnGroup2, MEDEXCEPTION);
1648   CPPUNIT_ASSERT_THROW(FIELD<double>::add(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
1649   CPPUNIT_ASSERT_THROW(FIELD<double>::sub(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
1650   CPPUNIT_ASSERT_THROW(FIELD<double>::addDeep(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
1651   CPPUNIT_ASSERT_THROW(FIELD<double>::subDeep(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
1652
1653   //CPPUNIT_ASSERT_NO_THROW();
1654   try
1655     {
1656       aPtr = FIELD<double>::mul(*aFieldOnGroup1, *aFieldOnGroup2);
1657       aPtr->removeReference();
1658       aPtr = FIELD<double>::div(*aFieldOnGroup1, *aFieldOnGroup2);
1659       aPtr->removeReference();
1660       aPtr = FIELD<double>::mulDeep(*aFieldOnGroup1, *aFieldOnGroup2);
1661       aPtr->removeReference();
1662       aPtr = FIELD<double>::divDeep(*aFieldOnGroup1, *aFieldOnGroup2);
1663       aPtr->removeReference();
1664
1665       *aFieldOnGroup1 *= *aFieldOnGroup2;
1666       *aFieldOnGroup1 /= *aFieldOnGroup2;
1667
1668       FIELD<double> *aPr = *aFieldOnGroup1 * *aFieldOnGroup2;
1669       FIELD<double> *aQu = *aFieldOnGroup1 / *aFieldOnGroup2;
1670       aPr->removeReference();
1671       aQu->removeReference();
1672     }
1673   catch (MEDEXCEPTION & ex)
1674     {
1675       CPPUNIT_FAIL(ex.what());
1676     }
1677   catch (...)
1678     {
1679       CPPUNIT_FAIL("Unknown exception");
1680     }
1681
1682   // restore MED units
1683   aFieldOnGroup1->setMEDComponentUnit(1, "unit1");
1684
1685   // check case of different operands: valueType
1686   //FIELD<int> * aFieldOnGroup4 =
1687   //  createFieldOnGroup<int>(aMeshOneMore, aGroupOneMore, "Test_Diff_Mesh", "test");
1688   //
1689   //CPPUNIT_ASSERT_THROW(FIELD<double>::add(*aFieldOnGroup4, *aFieldOnGroup3), MEDEXCEPTION);
1690   //CPPUNIT_ASSERT_THROW(*aFieldOnGroup4 - *aFieldOnGroup3, MEDEXCEPTION);
1691   //CPPUNIT_ASSERT_THROW(*aFieldOnGroup4 *= *aFieldOnGroup3, MEDEXCEPTION);
1692   //delete aFieldOnGroup4;
1693
1694   // check case of different operands: numberOfComponents
1695   //#ifdef ENABLE_FAULTS
1696   // (BUG) Cannot allocate value of higher dimension because of _componentsTypes reinitialization
1697   aFieldOnGroup1->deallocValue();
1698   //CPPUNIT_ASSERT_THROW(aFieldOnGroup1->allocValue(/*dim*/5), MEDEXCEPTION);
1699   CPPUNIT_ASSERT_NO_THROW(aFieldOnGroup1->allocValue(/*dim*/5));
1700   //#endif
1701   //#ifdef ENABLE_FORCED_FAILURES
1702   // YES THERE MUST BE AN EXCEPTION
1703   //   CPPUNIT_FAIL("Segmentation fault on attempt to allocate value of higher dimension."
1704   //                " Must be MEDEXCEPTION instead. And on attempt to change nb.components"
1705   //                " must be the same behaviour.");
1706   //#endif
1707   aFieldOnGroup1->setNumberOfComponents(5);
1708
1709   CPPUNIT_ASSERT_THROW(FIELD<double>::sub(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
1710   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 * *aFieldOnGroup2, MEDEXCEPTION);
1711   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 /= *aFieldOnGroup2, MEDEXCEPTION);
1712
1713   // check case of different operands: numberOfValues
1714   aFieldOnGroup1->deallocValue();
1715   aFieldOnGroup1->allocValue(2, nbVals + 1);
1716   // be carefull: aFieldOnGroup1 reallocated and contains random values
1717
1718   CPPUNIT_ASSERT_THROW(FIELD<double>::mul(*aFieldOnGroup1, *aFieldOnGroup2), MEDEXCEPTION);
1719   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 / *aFieldOnGroup2, MEDEXCEPTION);
1720   CPPUNIT_ASSERT_THROW(*aFieldOnGroup1 += *aFieldOnGroup2, MEDEXCEPTION);
1721
1722   // restore aFieldOnGroup1
1723   aFieldOnGroup1->deallocValue();
1724   aFieldOnGroup1->allocValue(2, nbVals);
1725   // be carefull: aFieldOnGroup1 reallocated and contains random values
1726
1727   aSubSupport1->removeReference();
1728   delete [] anElems1;
1729
1730   aScalarProduct->removeReference();
1731   aSubField1->removeReference();
1732   anAreaField->removeReference();
1733   barycenter->removeReference();
1734   aFieldOnGroup1->removeReference();
1735   aFieldOnGroup2->removeReference();
1736   aFieldOnGroup3->removeReference();
1737
1738   aMesh->removeReference();
1739   aMeshOneMore->removeReference();
1740
1741   /////////////////////
1742   // TEST 5: Drivers //
1743   /////////////////////
1744   testDrivers();
1745 }
1746
1747 /*!
1748  *  Check methods (2), defined in MEDMEM_FieldConvert.hxx:
1749  *  (+) template <class T> FIELD<T,FullInterlace> * FieldConvert(const FIELD<T,NoInterlace> & field);
1750  *  (+) template <class T> FIELD<T,NoInterlace> * FieldConvert(const FIELD<T,FullInterlace> & field);
1751  */
1752 void MEDMEMTest::testFieldConvert()
1753 {
1754   // create an empty integer field 2x10
1755   FIELD<int, FullInterlace> * aField_FING = new FIELD<int, FullInterlace>();
1756
1757   aField_FING->setName("Field_FING");
1758   aField_FING->setDescription("Field full interlace no gauss");
1759
1760   aField_FING->setNumberOfComponents(2);
1761   aField_FING->setNumberOfValues(10);
1762
1763   string aCompsNames[2] =
1764     {
1765       "Pos", "Neg"
1766     };
1767   string aCompsDescs[2] =
1768     {
1769       "+", "-"
1770     };
1771   string aMEDCompsUnits[2] =
1772     {
1773       "unit1", "unit2"
1774     };
1775   UNIT   aCompsUnits[2];
1776
1777   aCompsUnits[0] = UNIT("u1", "descr1");
1778   aCompsUnits[1] = UNIT("u2", "descr2");
1779
1780   aField_FING->setComponentsNames(aCompsNames);
1781   aField_FING->setComponentsDescriptions(aCompsDescs);
1782   aField_FING->setMEDComponentsUnits(aMEDCompsUnits);
1783   aField_FING->setComponentsUnits(aCompsUnits);
1784
1785   // check UNITs (for testField())
1786   const UNIT * aCompsUnitsBack = aField_FING->getComponentsUnits();
1787   CPPUNIT_ASSERT(aCompsUnits[0].getName() == aCompsUnitsBack[0].getName());
1788   CPPUNIT_ASSERT(aCompsUnits[1].getName() == aCompsUnitsBack[1].getName());
1789
1790   const UNIT * aCompUnitBack1 = aField_FING->getComponentUnit(1);
1791   const UNIT * aCompUnitBack2 = aField_FING->getComponentUnit(2);
1792   CPPUNIT_ASSERT(aCompsUnits[0].getName() == aCompUnitBack1->getName());
1793   CPPUNIT_ASSERT(aCompsUnits[1].getName() == aCompUnitBack2->getName());
1794
1795   // create one more field by copy
1796   FIELD<int, FullInterlace> * aField_FIGG = new FIELD<int, FullInterlace>(*aField_FING);
1797
1798   // values
1799   int values_FING[20] =
1800     {
1801       7,- 7, 14,-14, 21,-21, 28,-28, 35,-35,
1802       42,-42, 49,-49, 56,-56, 63,-63, 70,-70
1803     };
1804
1805   /////////////////////
1806   // TEST 1: NoGauss //
1807   /////////////////////
1808
1809   MEDMEM_ArrayInterface<int,FullInterlace,NoGauss>::Array * anArray_FING =
1810     new MEDMEM_ArrayInterface<int,FullInterlace,NoGauss>::Array
1811     (values_FING, /*dim*/2, /*nbelem*/10, /*shallowCopy*/false, /*ownershipOfValues*/false);
1812   aField_FING->setArray(anArray_FING);
1813   // no need to delete anArray_FING, because it will be deleted in destructor of aField_FING
1814
1815   // 1. FullInterlace -> NoInterlace
1816   FIELD<int, NoInterlace> * aField_NING = FieldConvert(*aField_FING);
1817   const int * values_NING = aField_NING->getValue();
1818
1819   for (int i = 0; i < 10; i++)
1820     {
1821       for (int j = 0; j < 2; j++)
1822         {
1823           CPPUNIT_ASSERT_EQUAL(values_FING[2*i + j], values_NING[10*j + i]);
1824         }
1825     }
1826
1827   // 2. NoInterlace -> FullInterlace
1828   FIELD<int, FullInterlace> * aField_FING_conv = FieldConvert(*aField_NING);
1829   const int * values_FING_conv = aField_FING_conv->getValue();
1830
1831   for (int i = 0; i < 10; i++)
1832     {
1833       for (int j = 0; j < 2; j++)
1834         {
1835           CPPUNIT_ASSERT_EQUAL(values_FING_conv[2*i + j], values_FING[2*i + j]);
1836           CPPUNIT_ASSERT_EQUAL(values_FING_conv[2*i + j], values_NING[10*j + i]);
1837         }
1838     }
1839
1840   aField_FING->removeReference();
1841   aField_NING->removeReference();
1842   aField_FING_conv->removeReference();
1843
1844   ///////////////////
1845   // TEST 2: Gauss //
1846   ///////////////////
1847   int nbelgeoc[2] =
1848     {
1849       1, 11
1850     };
1851   int nbgaussgeo[2] =
1852     {
1853       -1, 1
1854     };
1855   MEDMEM_ArrayInterface<int,FullInterlace,Gauss>::Array * anArray_FIGG =
1856     new MEDMEM_ArrayInterface<int,FullInterlace,Gauss>::Array
1857     (values_FING, /*dim*/2, /*nbelem*/10, /*nbtypegeo*/1, /*nbelgeoc*/nbelgeoc,
1858      /*nbgaussgeo*/nbgaussgeo, /*shallowCopy*/false, /*ownershipOfValues*/false);
1859   aField_FIGG->setArray(anArray_FIGG);
1860   // no need to delete anArray_FIGG, because it will be deleted in destructor of aField_FIGG
1861
1862   // 1. FullInterlace -> NoInterlace
1863   FIELD<int, NoInterlace> * aField_NIGG = FieldConvert(*aField_FIGG);
1864   const int * values_NIGG = aField_NIGG->getValue();
1865
1866   for (int i = 0; i < 10; i++)
1867     {
1868       for (int j = 0; j < 2; j++)
1869         {
1870           CPPUNIT_ASSERT_EQUAL(values_FING[2*i + j], values_NIGG[10*j + i]);
1871         }
1872     }
1873
1874   // 2. NoInterlace -> FullInterlace
1875   FIELD<int, FullInterlace> * aField_FIGG_conv = FieldConvert(*aField_NIGG);
1876   const int * values_FIGG_conv = aField_FIGG_conv->getValue();
1877
1878   for (int i = 0; i < 10; i++)
1879     {
1880       for (int j = 0; j < 2; j++)
1881         {
1882           CPPUNIT_ASSERT_EQUAL(values_FIGG_conv[2*i + j], values_FING[2*i + j]);
1883           CPPUNIT_ASSERT_EQUAL(values_FIGG_conv[2*i + j], values_NIGG[10*j + i]);
1884         }
1885     }
1886
1887   aField_FIGG->removeReference();
1888   aField_NIGG->removeReference();
1889   aField_FIGG_conv->removeReference();
1890   {
1891     // create an empty integer field 2x10
1892     FIELD<int, FullInterlace> * aField = new FIELD<int, FullInterlace>();
1893
1894     aField->setName("aField");
1895     aField->setDescription("Field full interlace no gauss");
1896
1897     aField->setNumberOfComponents(2);
1898     aField->setNumberOfValues(10);
1899
1900     aField->setComponentsNames(aCompsNames);
1901     aField->setComponentsDescriptions(aCompsDescs);
1902     aField->setMEDComponentsUnits(aMEDCompsUnits);
1903
1904     MEDMEM_ArrayInterface<int,FullInterlace,NoGauss>::Array * anArray =
1905       new MEDMEM_ArrayInterface<int,FullInterlace,NoGauss>::Array
1906       (values_FING, /*dim*/2, /*nbelem*/10, /*shallowCopy*/false, /*ownershipOfValues*/false);
1907     aField->setArray(anArray);
1908     // no need to delete anArray, because it will be deleted in destructor of aField
1909
1910     FIELD<int, NoInterlace> * aField_conv = FieldConvert(*aField);
1911     aField->removeReference();
1912     CPPUNIT_ASSERT(aField_conv);
1913     aField_conv->removeReference();
1914   }
1915 }
1916
1917 //================================================================================
1918 /*!
1919  * \brief 0020582: [CEA 368] MEDMEM don't work with a same field on NODES and CELLS
1920  * Before fixing the issue there was the error:
1921  *   RuntimeError: MED Exception in .../MED_SRC/src/MEDMEM/MEDMEM_MedFieldDriver22.hxx [503] : MED_FIELD_DRIVER<T>::createFieldSupportPart1(...) Field |field on NODEs and CELLs| with (ndt,or) = (-1,-1) must not be defined on different entity types
1922  */
1923 //================================================================================
1924
1925 void MEDMEMTest::testReadFieldOnNodesAndCells()
1926 {
1927   const string outfile = makeTmpFile("field_on_nodes_and_cells.med");
1928   const string fieldName = "field on NODEs and CELLs";
1929
1930   MEDMEMTest_TmpFilesRemover aTmpFilesRemover;
1931   aTmpFilesRemover.Register( outfile );
1932
1933   // write file with a field on NODEs and CELLs
1934
1935   using namespace MED_EN;
1936
1937   MEDMEM::MESH *mesh=MEDMEMTest_createTestMesh();
1938   int drv = mesh->addDriver( MED_DRIVER, outfile, mesh->getName());
1939   mesh->write(drv);
1940
1941   const SUPPORT *supportOnCells=mesh->getSupportOnAll( MED_CELL );
1942   const SUPPORT *supportOnNodes=mesh->getSupportOnAll( MED_NODE );
1943   supportOnCells->addReference();
1944   supportOnNodes->addReference();
1945   mesh->removeReference();
1946   int numberOfCells = supportOnCells->getNumberOfElements(MED_ALL_ELEMENTS);
1947   int numberOfNodes = supportOnNodes->getNumberOfElements(MED_ALL_ELEMENTS);
1948
1949   PointerOf<double> cellValues( numberOfCells ), nodeValues( numberOfNodes );
1950   for ( int i = 0; i < numberOfCells; ++i ) cellValues[i] = i;
1951   for ( int i = 0; i < numberOfNodes; ++i ) nodeValues[i] = -i;
1952
1953   FIELD<double> *wrFieldOnCells=new FIELD<double>(supportOnCells,1);
1954   wrFieldOnCells->setName(fieldName);
1955   wrFieldOnCells->setComponentName(1,"Vx");
1956   wrFieldOnCells->setComponentDescription(1,"comp1");
1957   wrFieldOnCells->setMEDComponentUnit(1,"unit1");
1958   wrFieldOnCells->setValue( cellValues );
1959   drv = wrFieldOnCells->addDriver(MED_DRIVER, outfile, fieldName);
1960   wrFieldOnCells->write( drv );
1961   wrFieldOnCells->removeReference();
1962
1963   FIELD<double> *wrFieldOnNodes=new FIELD<double>(supportOnNodes,1);
1964   wrFieldOnNodes->setName(fieldName);
1965   wrFieldOnNodes->setComponentName(1,"Vx");
1966   wrFieldOnNodes->setComponentDescription(1,"comp1");
1967   wrFieldOnNodes->setMEDComponentUnit(1,"unit1");
1968   wrFieldOnNodes->setValue( nodeValues );
1969   drv = wrFieldOnNodes->addDriver(MED_DRIVER, outfile, fieldName);
1970   wrFieldOnNodes->write( drv );
1971   wrFieldOnNodes->removeReference();
1972
1973   //  READ FIELDS BACK
1974
1975   //  field on CELLs
1976   FIELD<double> *cellField=new FIELD<double>(supportOnCells, MED_DRIVER, outfile, fieldName, -1, -1 );
1977   CPPUNIT_ASSERT_EQUAL( MED_CELL, cellField->getSupport()->getEntity());
1978   CPPUNIT_ASSERT_DOUBLES_EQUAL( numberOfCells-1, cellField->getValueIJ( numberOfCells, 1 ),1e-20);
1979   cellField->removeReference();
1980
1981   //  field on NODEs
1982   FIELD<double> *nodeField=new FIELD<double>(supportOnNodes, MED_DRIVER, outfile, fieldName, -1, -1 );
1983   CPPUNIT_ASSERT_EQUAL( MED_NODE, nodeField->getSupport()->getEntity());
1984   CPPUNIT_ASSERT_DOUBLES_EQUAL( -(numberOfNodes-1), nodeField->getValueIJ( numberOfNodes, 1 ),1e-20);
1985   nodeField->removeReference();
1986
1987   supportOnCells->removeReference();
1988   supportOnNodes->removeReference();
1989 }
1990
1991
1992 //================================================================================
1993 /*!
1994  * \brief 0021052: [CEA 431] Computation of Gauss point location
1995  *
1996  * Check method GetGaussPointsCoordinates() defined in MEDMEM_Field.hxx:
1997  */
1998 //================================================================================
1999 void MEDMEMTest::testGetGaussPointsCoordinates()
2000 {
2001   const int spaceDimension = 3;
2002   const int numberOfNodes = 56;
2003
2004   string names[3] =
2005     {
2006       "X","Y","Z"
2007     };
2008   string units[3] =
2009     {
2010       "cm","cm","cm"
2011     };
2012
2013   const int numberOfCellTypes = 14;
2014
2015   double coordinates [ spaceDimension*numberOfNodes ] =
2016     {
2017       0.0,     0.0,      0.0,   //N1
2018       0.5,     0.5,      0.5,   //N2
2019       1.0,     1.0,      1.0,   //N3
2020
2021       1.0,     1.0,      0.0,   //N4
2022       2.0,     2.5,      0.0,   //N5
2023       6.0,     1.5,      0.0,   //N6
2024       1.0,     2.0,      0.0,   //N7
2025       4.5,     2.5,      0.0,   //N8
2026       4.0,     0.5,      0.0,   //N9
2027
2028       0.0,     4.0,      0.0,   //N10
2029       4.0,     4.0,      0.0,   //N11
2030       4.0,     0.0,      0.0,   //N12
2031       0.0,     2.0,      0.0,   //N13
2032       2.0,     4.0,      0.0,   //N14
2033       4.0,     2.0,      0.0,   //N15
2034       2.0,     0.0,      0.0,   //N16
2035
2036       0.0,     6.0,      0.0,   //N17
2037       3.0,     3.0,      0.0,   //N18
2038       1.3,     3.0,      3.0,   //N19
2039
2040       0.0,     3.0,      0.0,   //N20
2041       1.5,     4.5,      0.0,   //N21
2042       1.5,     1.5,      0.0,   //N22
2043       0.65,    1.5,      1.5,   //N23
2044       0.65,    4.5,      1.5,   //N24
2045       2.15,    3.0,      1.5,   //N25
2046
2047       2.0,     2.0,      2.0,   //N26
2048       3.0,     1.0,      1.0,   //N27
2049       3.0,     3.0,      1.0,   //N28
2050       1.0,     3.0,      1.0,   //N29
2051       1.0,     1.0,      1.0,   //N30
2052
2053       0.0,     3.0,      0.0,   //N31
2054       2.0,     0.0,      0.0,   //N32
2055       0.0,     0.0,      6.0,   //N33
2056       0.0,     3.0,      6.0,   //N34
2057       3.0,     0.0,      6.0,   //N35
2058
2059       0.0,     1.5,      0.0,   //N36
2060       1.5,     1.5,      0.0,   //N37
2061       1.5,     0.0,      0.0,   //N38
2062       0.0,     1.5,      6.0,   //N39
2063       1.5,     1.5,      6.0,   //N40
2064       1.5,     0.0,      6.0,   //N41
2065       0.0,     0.0,      3.0,   //N42
2066       0.0,     3.0,      3.0,   //N43
2067       3.0,     0.0,      3.0,   //N44
2068
2069       0.0,     0.0,      4.0,   //N45
2070       0.0,     4.0,      4.0,   //N46
2071       4.0,     4.0,      4.0,   //N47
2072       4.0,     0.0,      4.0,   //N48
2073       0.0,     2.0,      4.0,   //N49
2074       2.0,     4.0,      4.0,   //N50
2075       4.0,     2.0,      4.0,   //N51
2076       2.0,     0.0,      4.0,   //N52
2077       0.0,     0.0,      2.0,   //N53
2078       0.0,     4.0,      2.0,   //N54
2079       4.0,     4.0,      2.0,   //N55
2080       4.0,     0.0,      2.0,   //N56
2081     };
2082
2083   MED_EN::medGeometryElement cellTypes[ numberOfCellTypes ] =
2084     {
2085       MED_EN::MED_SEG2,
2086       MED_EN::MED_SEG3,
2087       MED_EN::MED_TRIA3,
2088       MED_EN::MED_TRIA6,
2089       MED_EN::MED_QUAD4,
2090       MED_EN::MED_QUAD8,
2091       MED_EN::MED_TETRA4,
2092       MED_EN::MED_TETRA10,
2093       MED_EN::MED_PYRA5,
2094       MED_EN::MED_PYRA13,
2095       MED_EN::MED_PENTA6,
2096       MED_EN::MED_PENTA15,
2097       MED_EN::MED_HEXA8,
2098       MED_EN::MED_HEXA20
2099     };
2100
2101   const int numberOfCells[numberOfCellTypes] =
2102     {
2103       1, 1,                   //1D
2104       1, 1, 1, 1,             //2D
2105       1, 1, 1, 1, 1, 1, 1, 1  //3D
2106     };
2107
2108   //Seg2 Connectivity
2109   int seg2C [ 2 ] =
2110     {
2111       1,3
2112     };
2113
2114   //Seg3 Connectivity
2115   int seg3C [ 3 ] =
2116     {
2117       1,3,2
2118     };
2119
2120   //Tria3 Connectivity
2121   int tria3C [ 3 ] =
2122     {
2123       4,5,6
2124     };
2125
2126   //Tria6 Connectivity
2127   int tria6C [ 6 ] =
2128     {
2129       4,5,6,7,8,9
2130     };
2131
2132   //Quad4 Connectivity
2133   int quad4C [4] =
2134     {
2135       1, 10, 11, 12
2136     };
2137
2138   //Quad8 Connectivity
2139   int quad8C [8] =
2140     {
2141       1, 10, 11, 12, 13, 14, 15, 16
2142     };
2143
2144   //Tetra4 Connectivity
2145   int tetra4C [4] =
2146     {
2147       1, 17, 18, 19
2148     };
2149
2150   //Tetra10 Connectivity
2151   int tetra10C [10] =
2152     {
2153       1, 17, 18, 19, 20, 21, 22, 23, 24, 25
2154     };
2155
2156   //Pyra13 Connectivity
2157   int pyra5C [5] =
2158     {
2159       1, 12, 11, 10, 26
2160     };
2161
2162   //Pyra13 Connectivity
2163   int pyra13C [13] =
2164     {
2165       1, 12, 11, 10, 26, 16, 15, 14, 13, 27, 28, 29, 30
2166     };
2167
2168   //Penta6 Connectivity
2169   int penta6C [6] =
2170     {
2171       1, 31, 32, 33, 34, 35
2172     };
2173
2174   //Penta6 Connectivity
2175   int penta15C [15] =
2176     {
2177       1, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44
2178     };
2179
2180   //Hexa8 Connectivity
2181   int hexa8C[8] =
2182     {
2183       1, 10, 11, 12, 45, 46, 47, 48
2184     };
2185
2186   //Hexa20 Connectivity
2187   int hexa20C[20] =
2188     {
2189       1, 10, 11, 12, 45, 46, 47, 48, 13, 14, 15, 16, 49, 50, 51, 52, 53, 54, 55, 56
2190     };
2191
2192
2193
2194   MEDMEM::MESHING* mesh = new MEDMEM::MESHING;
2195   mesh->setCoordinates(spaceDimension, numberOfNodes, coordinates,
2196                        "CARTESIAN", MED_EN::MED_FULL_INTERLACE);
2197   mesh->setCoordinatesNames(names);
2198   mesh->setCoordinatesUnits(units);
2199
2200   //connectivities
2201   mesh->setNumberOfTypes(numberOfCellTypes, MED_EN::MED_CELL);
2202   mesh->setTypes(cellTypes, MED_EN::MED_CELL);
2203   mesh->setNumberOfElements(numberOfCells, MED_EN::MED_CELL);
2204
2205   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_SEG2, seg2C );
2206   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_SEG3, seg3C );
2207   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_TRIA3, tria3C );
2208   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_TRIA6, tria6C );
2209   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_QUAD4, quad4C );
2210   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_QUAD8, quad8C );
2211   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_TETRA4, tetra4C );
2212   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_TETRA10, tetra10C );
2213   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_PYRA5, pyra5C );
2214   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_PYRA13, pyra13C );
2215   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_PENTA6, penta6C );
2216   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_PENTA15, penta15C );
2217   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_HEXA8, hexa8C );
2218   mesh->setConnectivity( MED_EN::MED_CELL, MED_EN::MED_HEXA20, hexa20C );
2219
2220
2221   //Support definition
2222   const SUPPORT* sup = mesh->getSupportOnAll( MED_EN::MED_CELL );
2223
2224   //Create test field
2225   FIELD<int>* aField = new FIELD<int>(sup,1);
2226
2227   //Gauss Localization definition:
2228   double v[1] =
2229     {
2230       0.3
2231     };
2232   double v_2[2] =
2233     {
2234       0.3, 0.3
2235     };
2236   double v_3[3] =
2237     {
2238       0.3, 0.3, 0.3
2239     };
2240   double v_4[4] =
2241     {
2242       0.3, 0.3, 0.3, 0.3
2243     };
2244
2245   //              --------- Seg2 localization ---------
2246   //              Nb of the gauss points = 1
2247   double seg2CooRef[2] =
2248     {
2249       -1.0 , 1.0
2250     };
2251   double seg2CooGauss[1] =
2252     {
2253       0.2
2254     };
2255
2256   GAUSS_LOCALIZATION<>* aseg2L = new GAUSS_LOCALIZATION<>("Seg2 Gauss localization",
2257                                                           MED_EN::MED_SEG2,
2258                                                           1,
2259                                                           seg2CooRef,
2260                                                           seg2CooGauss,
2261                                                           v);
2262
2263   //              --------- Seg3 localization ---------
2264   //              Nb of the gauss points = 1
2265   double seg3CooRef[3] =
2266     {
2267       -1.0, 1.0, 0.0
2268     };
2269   double seg3CooGauss[1] =
2270     {
2271       0.2
2272     };
2273
2274   GAUSS_LOCALIZATION<>* aseg3L = new GAUSS_LOCALIZATION<>("Seg3 Gauss localization",
2275                                                           MED_EN::MED_SEG3,
2276                                                           1,
2277                                                           seg3CooRef,
2278                                                           seg3CooGauss,
2279                                                           v);
2280   //              --------- Tria3 localization ---------
2281   //              Nb of the gauss points = 2
2282   double tria3CooRef[6] =
2283     {
2284       0.0, 0.0, 1.0 , 0.0, 0.0, 1.0
2285     };
2286
2287   double tria3CooGauss[4] =
2288     {
2289       0.1, 0.8, 0.2, 0.7
2290     };
2291
2292   GAUSS_LOCALIZATION<>* atria3L = new GAUSS_LOCALIZATION<>("Tria3 Gauss localization",
2293                                                            MED_EN::MED_TRIA3,
2294                                                            2,
2295                                                            tria3CooRef,
2296                                                            tria3CooGauss,
2297                                                            v_2);
2298
2299   //              --------- Tria6 localization ---------
2300   //              Nb of the gauss points = 3
2301   double tria6CooRef[12] =
2302     {
2303       0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 0.5, 0.5, 0.0, 0.5
2304     };
2305
2306   double tria6CooGauss[6] =
2307     {
2308       0.3, 0.2, 0.2, 0.1, 0.2, 0.4
2309     };
2310
2311   GAUSS_LOCALIZATION<>* atria6L = new GAUSS_LOCALIZATION<>("Tria6 Gauss localization",
2312                                                            MED_EN::MED_TRIA6,
2313                                                            3,
2314                                                            tria6CooRef,
2315                                                            tria6CooGauss,
2316                                                            v_3);
2317   //              --------- Quad4 localization ---------
2318   //              Nb of the gauss points = 4
2319   double quad4CooRef[8] =
2320     {
2321       -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0, 1.0
2322     };
2323
2324   double quad4CooGauss[8] =
2325     { //Size is type/10*NbGauss = (204/100)*4 = 8
2326       0.3, 0.2, 0.2, 0.1, 0.2, 0.4, 0.15, 0.27
2327     };
2328
2329   GAUSS_LOCALIZATION<>* aquad4L = new GAUSS_LOCALIZATION<>("Quad8 Gauss localization",
2330                                                            MED_EN::MED_QUAD4,
2331                                                            4,
2332                                                            quad4CooRef,
2333                                                            quad4CooGauss,
2334                                                            v_4);
2335   //              --------- Quad8 localization ---------
2336   //              Nb of the gauss points = 4
2337   double quad8CooRef[16] =
2338     {
2339       -1.0, -1.0,
2340       1.0, -1.0,
2341       1.0,  1.0,
2342       -1.0,  1.0,
2343       0.0, -1.0,
2344       1.0,  0.0,
2345       0.0,  1.0,
2346       -1.0,  0.0
2347     };
2348
2349   double quad8CooGauss[8] =
2350     {
2351       0.34, 0.16, 0.21, 0.3, 0.23, 0.4, 0.14, 0.37
2352     };
2353
2354   GAUSS_LOCALIZATION<>* aquad8L = new GAUSS_LOCALIZATION<>("Quad8 Gauss localization",
2355                                                            MED_EN::MED_QUAD8,
2356                                                            4,
2357                                                            quad8CooRef,
2358                                                            quad8CooGauss,
2359                                                            v_4);
2360
2361   //              --------- Tetra4 localization
2362   //              Nb of the gauss points = 1
2363   double tetra4CooRef[12] =
2364     {
2365       0.0,  1.0,  0.0,
2366       0.0,  0.0,  1.0,
2367       0.0,  0.0,  0.0,
2368       1.0,  0.0,  0.0
2369     };
2370
2371   double tetra4CooGauss[3] =
2372     {
2373       0.34, 0.16, 0.21
2374     };
2375
2376   GAUSS_LOCALIZATION<>* atetra4L = new GAUSS_LOCALIZATION<>("Tetra4 Gauss localization",
2377                                                             MED_EN::MED_TETRA4,
2378                                                             1,
2379                                                             tetra4CooRef,
2380                                                             tetra4CooGauss,
2381                                                             v);
2382
2383   //              --------- Tetra10 localization
2384   //              Nb of the gauss points = 1
2385   double tetra10CooRef[30] =
2386     {
2387       0.0, 1.0, 0.0,
2388       0.0, 0.0, 0.0,
2389       0.0, 0.0, 1.0,
2390       1.0, 0.0, 0.0,
2391       0.0, 0.5, 0.0,
2392       0.0, 0.0, 0.5,
2393       0.0, 0.5, 0.5,
2394       0.5, 0.5, 0.0,
2395       0.5, 0.0, 0.0,
2396       0.5, 0.0, 0.5,
2397     };
2398
2399   double tetra10CooGauss[3] =
2400     {
2401       0.2, 0.3, 0.1
2402     };
2403
2404   GAUSS_LOCALIZATION<>* atetra10L = new GAUSS_LOCALIZATION<>("Tetra10 Gauss localization",
2405                                                              MED_EN::MED_TETRA10,
2406                                                              1,
2407                                                              tetra10CooRef,
2408                                                              tetra10CooGauss,
2409                                                              v);
2410   //              --------- Pyra5 localization
2411   //              Nb of the gauss points = 1
2412   double pyra5CooRef[15] =
2413     {
2414       1.0,  0.0,  0.0,
2415       0.0,  1.0,  0.0,
2416       -1.0,  0.0,  0.0,
2417       0.0, -1.0,  0.0,
2418       0.0,  0.0,  1.0
2419     };
2420
2421   double pyra5CooGauss[3] =
2422     {
2423       0.2, 0.3, 0.1
2424     };
2425
2426   GAUSS_LOCALIZATION<>* apyra5L = new GAUSS_LOCALIZATION<>("Pyra5 Gauss localization",
2427                                                            MED_EN::MED_PYRA5,
2428                                                            1,
2429                                                            pyra5CooRef,
2430                                                            pyra5CooGauss,
2431                                                            v);
2432
2433   //              --------- Pyra13 localization
2434   //              Nb of the gauss points = 1
2435   double pyra13CooRef[39] =
2436     {
2437       1.0,  0.0, 0.0,
2438       0.0,  1.0, 0.0,
2439       -1.0,  0.0, 0.0,
2440       0.0, -1.0, 0.0,
2441       0.0,  0.0, 1.0,
2442       0.5,  0.5, 0.0,
2443       -0.5,  0.5, 0.0,
2444       -0.5, -0.5, 0.0,
2445       0.5, -0.5, 0.0,
2446       0.5,  0.0, 0.5,
2447       0.0,  0.5, 0.5,
2448       -0.5,  0.0, 0.5,
2449       0.0, -0.5, 0.5
2450     };
2451
2452   double pyra13CooGauss[3] =
2453     {
2454       0.1, 0.2, 0.7
2455     };
2456
2457   GAUSS_LOCALIZATION<>* apyra13L = new GAUSS_LOCALIZATION<>("Pyra13 Gauss localization",
2458                                                             MED_EN::MED_PYRA13,
2459                                                             1,
2460                                                             pyra13CooRef,
2461                                                             pyra13CooGauss,
2462                                                             v);
2463   //              --------- Penta6 localization
2464   //              Nb of the gauss points = 1
2465   double penta6CooRef[18] =
2466     {
2467       -1.0,   1.0,   0.0,
2468       -1.0,  -0.0,   1.0,
2469       -1.0,   0.0,   0.0,
2470       1.0,   1.0,   0.0,
2471       1.0,   0.0,   1.0,
2472       1.0,   0.0,   0.0
2473     };
2474
2475   double penta6CooGauss[3] =
2476     {
2477       0.2, 0.3, 0.1
2478     };
2479
2480   GAUSS_LOCALIZATION<>* apenta6L = new GAUSS_LOCALIZATION<>("Penta6 Gauss localization",
2481                                                             MED_EN::MED_PENTA6,
2482                                                             1,
2483                                                             penta6CooRef,
2484                                                             penta6CooGauss,
2485                                                             v);
2486
2487   //              --------- Penta15 localization
2488   //              Nb of the gauss points = 1
2489   double penta15CooRef[45] =
2490     {
2491       -1.0,  1.0,   0.0,
2492       -1.0,  0.0,   1.0,
2493       -1.0,  0.0,   0.0,
2494       1.0,  1.0,   0.0,
2495       1.0,  0.0,   1.0,
2496       1.0,  0.0,   0.0,
2497       -1.0,  0.5,   0.5,
2498       -1.0,  0.0,   0.5,
2499       -1.0,  0.5,   0.0,
2500       0.0,  1.0,   0.0,
2501       0.0,  0.0,   1.0,
2502       0.0,  0.0,   0.0,
2503       1.0,  0.5,   0.5,
2504       1.0,  0.0,   0.5,
2505       1.0,  0.5,   0.0
2506     };
2507
2508   double penta15CooGauss[3] =
2509     {
2510       0.2, 0.3, 0.15
2511     };
2512
2513   GAUSS_LOCALIZATION<>* apenta15L = new GAUSS_LOCALIZATION<>("Penta15 Gauss localization",
2514                                                              MED_EN::MED_PENTA15,
2515                                                              1,
2516                                                              penta15CooRef,
2517                                                              penta15CooGauss,
2518                                                              v);
2519
2520   //              --------- Hexa8 localization
2521   //              Nb of the gauss points = 1
2522   double hexa8CooRef [24] =
2523     {
2524       -1.0,  -1.0,  -1.0,
2525       1.0,  -1.0,  -1.0,
2526       1.0,   1.0,  -1.0,
2527       -1.0,   1.0,  -1.0,
2528       -1.0,  -1.0,   1.0,
2529       1.0,  -1.0,   1.0,
2530       1.0,   1.0,   1.0,
2531       -1.0,   1.0,   1.0
2532     };
2533
2534   double hexa8CooGauss[3] =
2535     {
2536       0.2, 0.3, 0.15
2537     };
2538
2539   GAUSS_LOCALIZATION<>* ahexa8L = new GAUSS_LOCALIZATION<>("Hexa8 Gauss localization",
2540                                                            MED_EN::MED_HEXA8,
2541                                                            1,
2542                                                            hexa8CooRef,
2543                                                            hexa8CooGauss,
2544                                                            v);
2545
2546   //              --------- Hexa20 localization
2547   //              Nb of the gauss points = 1
2548   double hexa20CooRef[60] =
2549     {
2550       -1.0,  -1.0,  -1.0,
2551       1.0,  -1.0,  -1.0,
2552       1.0,   1.0,  -1.0,
2553       -1.0,   1.0,  -1.0,
2554       -1.0,  -1.0,   1.0,
2555       1.0,  -1.0,   1.0,
2556       1.0,   1.0,   1.0,
2557       -1.0,   1.0,   1.0,
2558       0.0,  -1.0,  -1.0,
2559       1.0,   0.0,  -1.0,
2560       0.0,   1.0,  -1.0,
2561       -1.0,   0.0,  -1.0,
2562       -1.0,  -1.0,   0.0,
2563       1.0,  -1.0,   0.0,
2564       1.0,   1.0,   0.0,
2565       -1.0,   1.0,   0.0,
2566       0.0,  -1.0,   1.0,
2567       1.0,   0.0,   1.0,
2568       0.0,   1.0,   1.0,
2569       -1.0,   0.0,   1.0
2570     };
2571
2572   double hexa20CooGauss[3] =
2573     {
2574       0.11, 0.3, 0.55
2575     };
2576
2577   GAUSS_LOCALIZATION<>* ahexa20L = new GAUSS_LOCALIZATION<>("Hexa20 Gauss localization",
2578                                                             MED_EN::MED_HEXA20,
2579                                                             1,
2580                                                             hexa20CooRef,
2581                                                             hexa20CooGauss,
2582                                                             v);
2583
2584
2585
2586   aField->setGaussLocalization(MED_EN::MED_SEG2, aseg2L);
2587   aField->setGaussLocalization(MED_EN::MED_SEG3, aseg3L);
2588   aField->setGaussLocalization(MED_EN::MED_TRIA3, atria3L);
2589   aField->setGaussLocalization(MED_EN::MED_TRIA6, atria6L);
2590   aField->setGaussLocalization(MED_EN::MED_QUAD4, aquad4L);
2591   aField->setGaussLocalization(MED_EN::MED_QUAD8, aquad8L);
2592   aField->setGaussLocalization(MED_EN::MED_TETRA4, atetra4L);
2593   aField->setGaussLocalization(MED_EN::MED_TETRA10, atetra10L);
2594   aField->setGaussLocalization(MED_EN::MED_PYRA5, apyra5L);
2595   aField->setGaussLocalization(MED_EN::MED_PYRA13, apyra13L);
2596   aField->setGaussLocalization(MED_EN::MED_PENTA6, apenta6L);
2597   aField->setGaussLocalization(MED_EN::MED_PENTA15, apenta15L);
2598   aField->setGaussLocalization(MED_EN::MED_HEXA8, ahexa8L);
2599   aField->setGaussLocalization(MED_EN::MED_HEXA20, ahexa20L);
2600
2601   FIELD<double>* aGaussCoords = aField->getGaussPointsCoordinates();
2602
2603
2604   //Coordinates to check result
2605   double seg2Coord[3] =
2606     {
2607       0.6,0.6,0.6
2608     };
2609   double seg3Coord[3] =
2610     {
2611       0.6,0.6,0.6
2612     };
2613
2614   double tria3Coord[3*2]  =
2615     {
2616       5.1, 1.55, 0.0,   //First GP Coordinates
2617       4.7, 1.65, 0.0    //Second GP Coordinates
2618     };
2619
2620   double tria6Coord[3*3] =
2621     {
2622       2.32, 1.52, 0.0, //First GP Coordinates
2623       1.6 , 1.32, 0.0, //Second GP Coordinates
2624       3.52, 1.26, 0.0  //Third GP Coordinates
2625     };
2626
2627   double quad4Coord[4*3] =
2628     {
2629       2.6, 1.6,  0.0,
2630       2.4, 1.8,  0.0,
2631       2.4, 1.2,  0.0,
2632       2.3, 1.46, 0.0
2633     };
2634
2635   double quad8Coord[4*3] =
2636     {
2637       2.32, 2.68,  0.0,
2638       2.6,  2.42,  0.0,
2639       2.8,  2.46,  0.0,
2640       2.74, 2.28,  0.0
2641     };
2642
2643   double tetra4Coord[3] =
2644     {
2645       1.312, 3.15, 1.02
2646     };
2647
2648   double tetra10Coord[3] =
2649     {
2650       0.56, 3.3, 0.6
2651     };
2652
2653   double pyra5Coord [3]=
2654     {
2655       2.18, 1.1, 0.2
2656     };
2657
2658   double pyra13Coord [3] =
2659     {
2660       1.18, 1.54, 0.98
2661     };
2662
2663   double penta6Coord [3] =
2664     {
2665       1.56, 0.3, 3.6
2666     };
2667
2668   double penta15Coord [3] =
2669     {
2670       1.613, 0.801, 4.374
2671     };
2672
2673   double hexa8Coord [3] =
2674     {
2675       2.6, 2.4, 2.3
2676     };
2677
2678   double hexa20Coord [3] =
2679     {
2680       2.31232, 2.3934, 1.55326
2681     };
2682
2683   //Check result of the calculation
2684   int ElemId = 1;
2685   double EPS = 0.000001;
2686   int idx = 0;
2687
2688   //Seg2:
2689   for(int j = 1; j<=3;j++)
2690     {
2691       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
2692                                    seg2Coord[j-1], EPS);
2693     }
2694
2695   //Seg3:
2696   ElemId++;
2697   for(int j = 1; j<=3;j++)
2698     {
2699       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
2700                                    seg3Coord[j-1], EPS);
2701     }
2702
2703   //Tria3
2704   ElemId++;
2705   for(int k = 1; k <=2 ; k++)
2706     {
2707       for(int j = 1; j<=3;j++)
2708         {
2709           CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,k),
2710                                        tria3Coord[idx++], EPS);
2711         }
2712     }
2713
2714   //Tria6
2715   ElemId++;
2716   idx=0;
2717   for(int k = 1; k <=3 ; k++)
2718     {
2719       for(int j = 1; j<=3;j++)
2720         {
2721           CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,k),
2722                                        tria6Coord[idx++], EPS);
2723         }
2724     }
2725
2726   //Quad4
2727   ElemId++;
2728   idx=0;
2729   for(int k = 1; k <=4 ; k++)
2730     {
2731       for(int j = 1; j<=3;j++)
2732         {
2733           CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,k),
2734                                        quad4Coord[idx++], EPS);
2735         }
2736     }
2737
2738   //Quad8
2739   ElemId++;
2740   idx=0;
2741   for(int k = 1; k <=4 ; k++)
2742     {
2743       for(int j = 1; j<=3;j++)
2744         {
2745           CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,k),
2746                                        quad8Coord[idx++], EPS);
2747         }
2748     }
2749
2750   //Tetra4
2751   ElemId++;
2752   for(int j = 1; j<=3;j++)
2753     {
2754       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
2755                                    tetra4Coord[j-1], EPS);
2756     }
2757
2758   //Tetra10
2759   ElemId++;
2760   for(int j = 1; j<=3;j++)
2761     {
2762       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
2763                                    tetra10Coord[j-1], EPS);
2764     }
2765
2766   //Pyra5
2767   ElemId++;
2768   for(int j = 1; j<=3;j++)
2769     {
2770       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
2771                                    pyra5Coord[j-1], EPS);
2772     }
2773
2774   //Penta15
2775   ElemId++;
2776   for(int j = 1; j<=3;j++)
2777     {
2778       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
2779                                    pyra13Coord[j-1], EPS);
2780     }
2781
2782   //Penta6
2783   ElemId++;
2784   for(int j = 1; j<=3;j++)
2785     {
2786       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
2787                                    penta6Coord[j-1], EPS);
2788     }
2789
2790   //Penta15
2791   ElemId++;
2792   for(int j = 1; j<=3;j++)
2793     {
2794       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
2795                                    penta15Coord[j-1], EPS);
2796     }
2797
2798   //Hexa8
2799   ElemId++;
2800   for(int j = 1; j<=3;j++)
2801     {
2802       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
2803                                    hexa8Coord[j-1], EPS);
2804     }
2805
2806   //Hexa20
2807   ElemId++;
2808   for(int j = 1; j<=3;j++)
2809     {
2810       CPPUNIT_ASSERT_DOUBLES_EQUAL(aGaussCoords->getValueIJK(ElemId,j,1),
2811                                    hexa20Coord[j-1], 0.0001);
2812     }
2813
2814   aGaussCoords->removeReference();
2815   aField->removeReference();
2816   mesh->removeReference();
2817 }