src/MEDMEMCppTest/MEDMEMTest_MeshAndMeshing_fault.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "MEDMEMTest.hxx"
  21 #include "MEDMEM_Meshing.hxx"
  22 #include "MEDMEM_Group.hxx"
  23 #include "MEDMEM_define.hxx"
  24 #include "MEDMEM_MedMeshDriver.hxx"
  25 #include "MEDMEM_Field.hxx"
  26 #include "MEDMEM_Grid.hxx"
  27
  28 #include <cppunit/Message.h>
  29 #include <cppunit/TestAssert.h>
  30
  31 #include <sstream>
  32 #include <cmath>
  33
  34 // use this define to enable lines, execution of which leads to Segmentation Fault
  35 //#define ENABLE_FAULTS
  36
  37 // use this define to enable CPPUNIT asserts and fails, showing bugs
  38 //#define ENABLE_FORCED_FAILURES
  39
  40 using namespace std;
  41 using namespace MEDMEM;
  42 using namespace MED_EN;
  43
  44 double dmax(double x, double y)
  45 {
  46 return (x>y)?x:y;
  47 }
  48
  49 double dmin(double x, double y)
  50 {
  51 return (x>y)?y:x;
  52 }
  53
  54 /*!
  55  *  Check methods (18), defined in MEDMEM_Meshing.hxx:
  56  *  class MESHING: public MESH
  57 {
  58  *   (+) MESHING();
  59  *   (+) ~MESHING();
  60  *   (+) void setSpaceDimension (const int SpaceDimension);
  61  *   (+) void setMeshDimension (const int MeshDimension);
  62  *   (+) void setNumberOfNodes (const int NumberOfNodes);
  63  *   (+) void setCoordinates (const int SpaceDimension, const int NumberOfNodes,
  64  *                                const double * Coordinates,
  65  *                                const string System, const MED_EN::medModeSwitch Mode);
  66  *   (+) void setCoordinatesSystem(const string System) throw (MEDEXCEPTION);
  67  *   (+) void setCoordinatesNames (const string * names);
  68  *   (+) void setCoordinateName (const string name, const int i);
  69  *   (+) void setCoordinatesUnits (const string * units);
  70  *   (+) void setCoordinateUnit (const string unit, const int i);
  71  *   (+) void setNumberOfTypes (const int NumberOfTypes,
  72  *                                  const MED_EN::medEntityMesh Entity) throw (MEDEXCEPTION);
  73  *   (+) void setTypes (const MED_EN::medGeometryElement * Types,
  74  *                          const MED_EN::medEntityMesh Entity) throw (MEDEXCEPTION);
  75  *   (+) void setNumberOfElements (const int * NumberOfElements,
  76  *                                 const MED_EN::medEntityMesh Entity) throw (MEDEXCEPTION);
  77  *   (+) void setConnectivity (const int * Connectivity, const MED_EN::medEntityMesh Entity,
  78  *                             const MED_EN::medGeometryElement Type) throw (MEDEXCEPTION);
  79  *   (+) void setPolygonsConnectivity (const int * ConnectivityIndex, const int * ConnectivityValue,
  80  *                                     int nbOfPolygons,
  81  *                                     const MED_EN::medEntityMesh Entity) throw (MEDEXCEPTION);
  82  *   (+) void setPolyhedraConnectivity (const int * PolyhedronIndex, const int * FacesIndex,
  83  *                                      const int * Nodes, int nbOfPolyhedra,
  84  *                                      const MED_EN::medEntityMesh Entity) throw (MEDEXCEPTION);
  85  *   (NOT YET IMPLEMENTED!!!) void setConnectivities (const int * ConnectivityIndex,
  86  *                                   const int * ConnectivityValue,
  87  *                                   const MED_EN::medConnectivity ConnectivityType,
  88  *                                   const MED_EN::medEntityMesh Entity) throw (MEDEXCEPTION);
  89  *   (+) void addGroup (const GROUP & Group) throw (MEDEXCEPTION);
  90  *
  91 }
  92  */
  93
  94 /*!
  95  *  Check methods (87), defined in MEDMEM_Mesh.hxx:
  96  *  class MESH : public RCBASE
  97 {
  98  *   (+) void init();
  99  *   (+) MESH();
 100  *   (+) MESH(MESH &m);
 101  *   (+) MESH & operator=(const MESH &m);
 102  *   (+) virtual bool operator==(const MESH& other) const;
 103  *   (+) virtual bool deepCompare(const MESH& other) const;
 104  *   (+) MESH(driverTypes driverType, const string & fileName="",
 105  *                const string & meshName="") throw (MEDEXCEPTION);
 106  *   (+) virtual ~MESH();
 107  *   (+) friend ostream & operator<<(ostream &os, const MESH &my);
 108  *   (+) int  addDriver(driverTypes driverType,
 109  *                          const string & fileName="Default File Name.med",
 110  *                          const string & driverName="Default Mesh Name",
 111  *                          MED_EN::med_mode_acces access=MED_EN::MED_REMP);
 112  *   (+) int  addDriver(GENDRIVER & driver);
 113  *   (+) void rmDriver(int index=0);
 114  *   (+) virtual void read(int index=0);
 115  *   (must be private) inline void read(const GENDRIVER & genDriver);
 116  *   (+) inline void write(int index=0, const string & driverName = "");
 117  *   (must be private) inline void write(const GENDRIVER & genDriver);
 118  *   (+) inline void setName(string name);
 119  *   (+) inline void setDescription(string description);
 120  *   (+) inline string getName() const;
 121  *   (+) inline string getDescription() const;
 122  *   (+) inline int getSpaceDimension() const;
 123  *   (+) inline int getMeshDimension() const;
 124  *   (+) inline bool getIsAGrid();
 125  *   (+) inline int getNumberOfNodes() const;
 126  *   (+) virtual inline const COORDINATE * getCoordinateptr() const;
 127  *   (+) inline string                     getCoordinatesSystem() const;
 128  *   (+) virtual inline const double *     getCoordinates(MED_EN::medModeSwitch Mode) const;
 129  *   (+) virtual inline const double       getCoordinate(int Number,int Axis) const;
 130  *   (+) inline const string *             getCoordinatesNames() const;
 131  *   (+) inline const string *             getCoordinatesUnits() const;
 132  *   (+) virtual inline int getNumberOfTypes(MED_EN::medEntityMesh Entity) const;
 133  *   (+) virtual int getNumberOfTypesWithPoly(MED_EN::medEntityMesh Entity) const;
 134  *   (+) virtual inline const MED_EN::medGeometryElement * getTypes(MED_EN::medEntityMesh Entity) const;
 135  *   (+) virtual MED_EN::medGeometryElement * getTypesWithPoly(MED_EN::medEntityMesh Entity) const;
 136  *   (+) virtual inline const CELLMODEL * getCellsTypes(MED_EN::medEntityMesh Entity) const;
 137  *   (+) virtual inline string * getCellTypeNames(MED_EN::medEntityMesh Entity) const;
 138  *   (+) virtual const int * getGlobalNumberingIndex(MED_EN::medEntityMesh Entity) const;
 139  *   (+) virtual inline int getNumberOfElements(MED_EN::medEntityMesh Entity,
 140  *                                        MED_EN::medGeometryElement Type) const;
 141  *   (+) virtual int getNumberOfElementsWithPoly(MED_EN::medEntityMesh Entity,
 142  *                                        MED_EN::medGeometryElement Type) const;
 143  *   (+) virtual inline bool existConnectivity(MED_EN::medConnectivity ConnectivityType,
 144  *                                       MED_EN::medEntityMesh Entity) const;
 145  *   (+) inline bool existPolygonsConnectivity(MED_EN::medConnectivity ConnectivityType,
 146  *                                       MED_EN::medEntityMesh Entity) const;
 147  *   (+) inline bool existPolyhedronConnectivity(MED_EN::medConnectivity ConnectivityType,
 148  *                                         MED_EN::medEntityMesh Entity) const;
 149  *   (+) virtual inline MED_EN::medGeometryElement getElementType
 150  *           (MED_EN::medEntityMesh Entity, int Number) const;
 151  *   (+) virtual inline MED_EN::medGeometryElement getElementTypeWithPoly
 152  *           (MED_EN::medEntityMesh Entity, int Number) const;
 153  *   (+) virtual inline void calculateConnectivity(MED_EN::medModeSwitch Mode,
 154  *                                            MED_EN::medConnectivity ConnectivityType,
 155  *                                            MED_EN::medEntityMesh Entity) const;
 156  *   (+) virtual inline int getConnectivityLength(MED_EN::medModeSwitch Mode,
 157  *                                             MED_EN::medConnectivity ConnectivityType,
 158  *                                             MED_EN::medEntityMesh Entity,
 159  *                                             MED_EN::medGeometryElement Type) const;
 160  *   (+) virtual inline const int * getConnectivity(MED_EN::medModeSwitch Mode,
 161  *                                             MED_EN::medConnectivity ConnectivityType,
 162  *                                             MED_EN::medEntityMesh Entity,
 163  *                                             MED_EN::medGeometryElement Type) const;
 164  *   (+) virtual inline const int * getConnectivityIndex
 165  *           (MED_EN::medConnectivity ConnectivityType, MED_EN::medEntityMesh Entity) const;
 166  *   (+) inline int getPolygonsConnectivityLength
 167  *           (MED_EN::medConnectivity ConnectivityType, MED_EN::medEntityMesh Entity) const;
 168  *   (+) inline const int * getPolygonsConnectivity
 169  *           (MED_EN::medConnectivity ConnectivityType, MED_EN::medEntityMesh Entity) const;
 170  *   (+) inline const int * getPolygonsConnectivityIndex
 171  *           (MED_EN::medConnectivity ConnectivityType, MED_EN::medEntityMesh Entity) const;
 172  *   (+) inline int getNumberOfPolygons(MED_EN::medEntityMesh Entity=MED_EN::MED_ALL_ENTITIES) const;
 173  *   (+) inline int getPolyhedronConnectivityLength(MED_EN::medConnectivity ConnectivityType) const;
 174  *   (+) inline const int * getPolyhedronConnectivity(MED_EN::medConnectivity ConnectivityType) const;
 175  *   (+) inline const int * getPolyhedronFacesIndex() const;
 176  *   (+) inline const int * getPolyhedronIndex(MED_EN::medConnectivity ConnectivityType) const;
 177  *   (+) inline int getNumberOfPolyhedronFaces() const;
 178  *   (+) inline int getNumberOfPolyhedron() const;
 179  *   (+) virtual int getElementNumber(MED_EN::medConnectivity ConnectivityType,
 180  *                                    MED_EN::medEntityMesh Entity, MED_EN::medGeometryElement Type,
 181  *                                    int * connectivity) const;
 182  *   (+) virtual inline int getReverseConnectivityLength
 183  *           (MED_EN::medConnectivity ConnectivityType, MED_EN::medEntityMesh Entity=MED_EN::MED_CELL) const;
 184  *   (+) virtual inline const int * getReverseConnectivity
 185  *           (MED_EN::medConnectivity ConnectivityType, MED_EN::medEntityMesh Entity=MED_EN::MED_CELL) const;
 186  *   (+) virtual inline int getReverseConnectivityIndexLength
 187  *           (MED_EN::medConnectivity ConnectivityType, MED_EN::medEntityMesh Entity=MED_EN::MED_CELL) const;
 188  *   (+) virtual inline const int * getReverseConnectivityIndex
 189  *           (MED_EN::medConnectivity ConnectivityType, MED_EN::medEntityMesh Entity=MED_EN::MED_CELL) const;
 190  *   (+) virtual int getNumberOfFamilies(MED_EN::medEntityMesh Entity) const;
 191  *   (+) virtual inline const vector<FAMILY*> getFamilies(MED_EN::medEntityMesh Entity) const;
 192  *   (+) virtual inline const FAMILY* getFamily(MED_EN::medEntityMesh Entity,int i) const;
 193  *   (+) virtual int getNumberOfGroups(MED_EN::medEntityMesh Entity) const;
 194  *   (+) virtual inline const vector<GROUP*> getGroups(MED_EN::medEntityMesh Entity) const;
 195  *   (+) virtual inline const GROUP* getGroup(MED_EN::medEntityMesh Entity,int i) const;
 196  *   (+) virtual inline const CONNECTIVITY* getConnectivityptr() const;
 197  *   (+) virtual SUPPORT * getBoundaryElements(MED_EN::medEntityMesh Entity) throw (MEDEXCEPTION);
 198  *   (+) SUPPORT * getSkin(const SUPPORT * Support3D) throw (MEDEXCEPTION);
 199  *   (+) virtual FIELD<double>* getVolume (const SUPPORT * Support) const throw (MEDEXCEPTION);
 200  *   (+) virtual FIELD<double>* getArea (const SUPPORT * Support) const throw (MEDEXCEPTION);
 201  *   (+) virtual FIELD<double>* getLength (const SUPPORT * Support) const throw (MEDEXCEPTION);
 202  *   (+) virtual FIELD<double>* getNormal (const SUPPORT * Support) const throw (MEDEXCEPTION);
 203  *   (+) virtual FIELD<double>* getBarycenter (const SUPPORT * Support) const throw (MEDEXCEPTION);
 204  *   (+) static SUPPORT * mergeSupports(const vector<SUPPORT *> Supports) throw (MEDEXCEPTION);
 205  *   (+) static SUPPORT * intersectSupports(const vector<SUPPORT *> Supports) throw (MEDEXCEPTION);
 206  *   (+) void createFamilies();
 207  *   (+) SUPPORT *buildSupportOnNodeFromElementList
 208  *            (const list<int>& listOfElt, MED_EN::medEntityMesh entity) const throw (MEDEXCEPTION);
 209  *   (do the same thing as buildSupportOnNodeFromElementList except that a SUPPORT is not created) void fillSupportOnNodeFromElementList
 210  *            (const list<int>& listOfElt, SUPPORT *supportToFill) const throw (MEDEXCEPTION);
 211  *   (+) SUPPORT *buildSupportOnElementsFromElementList
 212  *            (const list<int>& listOfElt, MED_EN::medEntityMesh entity) const throw (MEDEXCEPTION);
 213  *   (+) int getElementContainingPoint(const double *coord);
 214  *   (+) vector< vector<double> > getBoundingBox() const;
 215  *   (+) template<class T> static
 216  *       FIELD<T> * mergeFields(const vector< FIELD<T> * > & others, bool meshCompare=false);
 217  *   (Presently disconnected in C++) virtual void addReference() const;
 218  *   (Presently disconnected in C++) virtual void removeReference() const;
 219  *
 220 }
 221  */
 222 void MEDMEMTest_testMeshAndMeshing()
 223 {
 224   string filename = getResourceFile("pointe.med");
 225   string meshname = "maa1";
 226   string filenameout21 = makeTmpFile("myMeshWrite4_pointe21.med");
 227
 228   // To remove tmp files from disk
 229   MEDMEMTest_TmpFilesRemover aRemover;
 230   aRemover.Register(filenameout21);
 231
 232   ////////////
 233   // TEST 1 //
 234   ////////////
 235
 236   MESH * myMesh= new MESH();
 237   myMesh->setName("FIRST_MESH");
 238   CPPUNIT_ASSERT(myMesh != NULL);
 239
 240   //test operator <<
 241   //#ifdef ENABLE_FAULTS
 242   //CPPUNIT_ASSERT_NO_THROW(cout << *myMesh << endl);
 243   //#endif
 244   //#ifdef ENABLE_FORCED_FAILURES
 245   //  CPPUNIT_FAIL("ERROR: operator << : if mesh is empty then attempt"
 246   //               " to get values from null object causes error");
 247   //#endif
 248
 249   //test operator =
 250   MESH myMesh1 = *myMesh;
 251
 252   //deepCompare
 253   //#ifdef ENABLE_FAULTS
 254   bool isEqual = false;
 255   CPPUNIT_ASSERT_NO_THROW(isEqual = myMesh1.deepCompare(*myMesh));
 256   CPPUNIT_ASSERT(isEqual);
 257   //#endif
 258   //#ifdef ENABLE_FORCED_FAILURES
 259   //  CPPUNIT_FAIL("ERROR: deepCompare(...) fails if mesh is empty");
 260   //#endif
 261
 262   //ensure it imposible to compare meshes
 263   MESH *myMeshPointer =  myMesh;
 264   //test operator ==
 265   CPPUNIT_ASSERT(*myMeshPointer == *myMesh);
 266
 267   delete myMesh;
 268
 269   //set a MESH object
 270   MESHING *myMeshing=new MESHING;
 271   myMeshing->setName("meshing");
 272   // define coordinates
 273
 274   int SpaceDimension = 3;
 275   int NumberOfNodes = 19;
 276   double Coordinates[57] =
 277     {
 278       0.0, 0.0, 0.0,
 279       0.0, 0.0, 1.0,
 280       2.0, 0.0, 1.0,
 281       0.0, 2.0, 1.0,
 282       -2.0, 0.0, 1.0,
 283       0.0, -2.0, 1.0,
 284       1.0, 1.0, 2.0,
 285       -1.0, 1.0, 2.0,
 286       -1.0, -1.0, 2.0,
 287       1.0, -1.0, 2.0,
 288       1.0, 1.0, 3.0,
 289       -1.0, 1.0, 3.0,
 290       -1.0, -1.0, 3.0,
 291       1.0, -1.0, 3.0,
 292       1.0, 1.0, 4.0,
 293       -1.0, 1.0, 4.0,
 294       -1.0, -1.0, 4.0,
 295       1.0, -1.0, 4.0,
 296       0.0, 0.0, 5.0
 297     };
 298   try
 299     {
 300       myMeshing->setCoordinates(SpaceDimension,NumberOfNodes,Coordinates,"CARTESIAN",MED_FULL_INTERLACE);
 301     }
 302   catch (const std::exception &e)
 303     {
 304       CPPUNIT_FAIL(e.what());
 305     }
 306   catch (...)
 307     {
 308       CPPUNIT_FAIL("Unknown exception");
 309     }
 310
 311   string Names[3] =
 312     {
 313       "X","Y","Z"
 314     };
 315   try
 316     {
 317       myMeshing->setCoordinatesNames(Names);
 318     }
 319   catch (const std::exception &e)
 320     {
 321       CPPUNIT_FAIL(e.what());
 322     }
 323   catch (...)
 324     {
 325       CPPUNIT_FAIL("Unknown exception");
 326     }
 327
 328   string Units[3] =
 329     {
 330       "cm","cm","cm"
 331     };
 332   try
 333     {
 334       myMeshing->setCoordinatesUnits(Units);
 335     }
 336   catch (const std::exception &e)
 337     {
 338       CPPUNIT_FAIL(e.what());
 339     }
 340   catch (...)
 341     {
 342       CPPUNIT_FAIL("Unknown exception");
 343     }
 344
 345   // define conectivities
 346
 347   // cell part
 348
 349   const int NumberOfTypes = 3;
 350   medGeometryElement Types[NumberOfTypes] =
 351     {
 352       MED_TETRA4,MED_PYRA5,MED_HEXA8
 353     };
 354   const int NumberOfElements[NumberOfTypes] =
 355     {
 356       12,2,2
 357     };
 358
 359   CPPUNIT_ASSERT_NO_THROW(myMeshing->setNumberOfTypes(NumberOfTypes,MED_CELL));
 360
 361   CPPUNIT_ASSERT_NO_THROW(myMeshing->setTypes(Types,MED_CELL));
 362
 363   CPPUNIT_ASSERT_NO_THROW(myMeshing->setNumberOfElements(NumberOfElements,MED_CELL));
 364
 365   const int sizeTetra = 12*4;
 366   int ConnectivityTetra[sizeTetra]=
 367     {
 368       1,2,3,6,
 369       1,2,4,3,
 370       1,2,5,4,
 371       1,2,6,5,
 372       2,7,4,3,
 373       2,8,5,4,
 374       2,9,6,5,
 375       2,10,3,6,
 376       2,7,3,10,
 377       2,8,4,7,
 378       2,9,5,8,
 379       2,10,6,9
 380     };
 381
 382   CPPUNIT_ASSERT_NO_THROW(myMeshing->setConnectivity( MED_CELL, MED_TETRA4, ConnectivityTetra ));
 383
 384   int ConnectivityPyra[2*5]=
 385     {
 386       7,8,9,10,2,
 387       15,18,17,16,19
 388     };
 389
 390   CPPUNIT_ASSERT_NO_THROW(myMeshing->setConnectivity( MED_CELL, MED_PYRA5, ConnectivityPyra ));
 391
 392   int ConnectivityHexa[2*8]=
 393     {
 394       11,12,13,14,7,8,9,10,
 395       15,16,17,18,11,12,13,14
 396     };
 397
 398   CPPUNIT_ASSERT_NO_THROW(myMeshing->setConnectivity( MED_CELL, MED_HEXA8, ConnectivityHexa ));
 399
 400   // face part
 401   const int NumberOfFacesTypes = 2;
 402   medGeometryElement FacesTypes[NumberOfFacesTypes] =
 403     {
 404       MED_TRIA3,MED_QUAD4
 405     };
 406   const int NumberOfFacesElements[NumberOfFacesTypes] =
 407     {
 408       4,4
 409     };
 410
 411   CPPUNIT_ASSERT_NO_THROW(myMeshing->setNumberOfTypes(NumberOfFacesTypes,MED_FACE));
 412   CPPUNIT_ASSERT_NO_THROW(myMeshing->setTypes(FacesTypes,MED_FACE));
 413   CPPUNIT_ASSERT_NO_THROW(myMeshing->setNumberOfElements(NumberOfFacesElements,MED_FACE));
 414   const int nbTria = 4;
 415   const int sizeTria = nbTria*3;
 416   int ConnectivityTria[sizeTria]=
 417     {
 418       1,4,3,
 419       1,5,4,
 420       1,6,5,
 421       1,3,6
 422     };
 423
 424   CPPUNIT_ASSERT_NO_THROW(myMeshing->setConnectivity( MED_FACE, MED_TRIA3, ConnectivityTria ));
 425   const int nbQua = 4;
 426   int ConnectivityQua[nbQua*4]=
 427     {
 428       7,8,9,10,
 429       11,12,13,14,
 430       11,7,8,12,
 431       12,8,9,13
 432     };
 433
 434   CPPUNIT_ASSERT_NO_THROW(myMeshing->setConnectivity( MED_FACE, MED_QUAD4, ConnectivityQua ));
 435
 436   int meshDimension = SpaceDimension; // because there 3D cells in the mesh
 437
 438   // edge part
 439
 440   // not yet implemented : if set, results are unpredictable.
 441
 442   // Some groups :
 443
 444   // Node :
 445   {
 446     GROUP *myGroup=new GROUP;
 447     myGroup->setName("SomeNodes");
 448     myGroup->setMesh(myMeshing);
 449     myGroup->setEntity(MED_NODE);
 450     myGroup->setNumberOfGeometricType(1);
 451     medGeometryElement myTypes[1] =
 452       {
 453         MED_NONE
 454       };
 455     myGroup->setGeometricType(myTypes);
 456     const int myNumberOfElements[1] =
 457       {
 458         4
 459       };
 460     myGroup->setNumberOfElements(myNumberOfElements);
 461     const int index[1+1] =
 462       {
 463         1,5
 464       };
 465     const int value[4]=
 466       {
 467         1,4,5,7
 468       };
 469     myGroup->setNumber(index,value);
 470     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
 471     myGroup->removeReference();
 472   }
 473   {
 474     GROUP *myGroup=new GROUP;
 475     myGroup->setName("OtherNodes");
 476     myGroup->setMesh(myMeshing);
 477     myGroup->setEntity(MED_NODE);
 478     myGroup->setNumberOfGeometricType(1);
 479     medGeometryElement myTypes[1] =
 480       {
 481         MED_NONE
 482       };
 483     myGroup->setGeometricType(myTypes);
 484     const int myNumberOfElements[1] =
 485       {
 486         3
 487       };
 488     myGroup->setNumberOfElements(myNumberOfElements);
 489     const int index[1+1] =
 490       {
 491         1,4
 492       };
 493     const int value[3]=
 494       {
 495         2,3,6
 496       };
 497     myGroup->setNumber(index,value);
 498     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
 499     myGroup->removeReference();
 500   }
 501
 502   // Cell :
 503   {
 504     GROUP *myGroup=new GROUP;
 505     myGroup->setName("SomeCells");
 506     myGroup->setMesh(myMeshing);
 507     myGroup->setEntity(MED_CELL);
 508     myGroup->setNumberOfGeometricType(3);
 509     medGeometryElement myTypes[3] =
 510       {
 511         MED_TETRA4,MED_PYRA5,MED_HEXA8
 512       };
 513     myGroup->setGeometricType(myTypes);
 514     const int myNumberOfElements[3] =
 515       {
 516         4,1,2
 517       };
 518     myGroup->setNumberOfElements(myNumberOfElements);
 519     const int index[3+1] =
 520       {
 521         1,5,6,8
 522       };
 523     const int value[4+1+2]=
 524       {
 525         2,7,8,12,
 526         13,
 527         15,16
 528       };
 529     myGroup->setNumber(index,value);
 530     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
 531     myGroup->removeReference();
 532   }
 533   {
 534     GROUP *myGroup=new GROUP;
 535     myGroup->setName("OtherCells");
 536     myGroup->setMesh(myMeshing);
 537     myGroup->setEntity(MED_CELL);
 538     myGroup->setNumberOfGeometricType(2);
 539     medGeometryElement myTypes[] =
 540       {
 541         MED_TETRA4,MED_PYRA5
 542       };
 543     myGroup->setGeometricType(myTypes);
 544     const int myNumberOfElements[] =
 545       {
 546         4,1
 547       };
 548     myGroup->setNumberOfElements(myNumberOfElements);
 549     const int index[2+1] =
 550       {
 551         1,5,6
 552       };
 553     const int value[4+1]=
 554       {
 555         3,4,5,9,
 556         14
 557       };
 558     myGroup->setNumber(index,value);
 559     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
 560     myGroup->removeReference();
 561   }
 562
 563   // Face :
 564   {
 565     GROUP *myGroup=new GROUP;
 566     myGroup->setName("SomeFaces");
 567     myGroup->setMesh(myMeshing);
 568     myGroup->setEntity(MED_FACE);
 569     myGroup->setNumberOfGeometricType(2);
 570     medGeometryElement myTypes[2] =
 571       {
 572         MED_TRIA3,MED_QUAD4
 573       };
 574     myGroup->setGeometricType(myTypes);
 575     const int myNumberOfElements[2] =
 576       {
 577         2,3
 578       };
 579     myGroup->setNumberOfElements(myNumberOfElements);
 580     const int index[2+1] =
 581       {
 582         1,3,6
 583       };
 584     const int value[2+3]=
 585       {
 586         2,4,
 587         5,6,8
 588       };
 589     myGroup->setNumber(index,value);
 590     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
 591     myGroup->removeReference();
 592   }
 593   {
 594     GROUP *myGroup=new GROUP;
 595     myGroup->setName("OtherFaces");
 596     myGroup->setMesh(myMeshing);
 597     myGroup->setEntity(MED_FACE);
 598     myGroup->setNumberOfGeometricType(1);
 599     medGeometryElement myTypes[1] =
 600       {
 601         MED_TRIA3
 602       };
 603     myGroup->setGeometricType(myTypes);
 604     const int myNumberOfElements[1] =
 605       {
 606         2
 607       };
 608     myGroup->setNumberOfElements(myNumberOfElements);
 609     const int index[1+1] =
 610       {
 611         1,3
 612       };
 613     const int value[2]=
 614       {
 615         1,3
 616       };
 617     myGroup->setNumber(index,value);
 618     CPPUNIT_ASSERT_NO_THROW(myMeshing->addGroup(*myGroup));
 619     myGroup->removeReference();
 620   }
 621
 622   //test Mesh(MESH &m)
 623   {
 624     MESH * myMesh2 = new MESH( *myMeshing );
 625     CPPUNIT_ASSERT(myMesh2->deepCompare(*myMeshing));
 626
 627     cout<<*myMesh2<<endl;
 628     ostringstream os;
 629     os << * myMesh2;
 630     CPPUNIT_ASSERT(os.str() != "");
 631
 632     CPPUNIT_ASSERT_EQUAL(myMesh2->getName(),(string)"meshing");
 633     CPPUNIT_ASSERT((myMesh2->getDescription()).size() == 0);
 634     myMesh2->setDescription("This class contains all information related to a 'meshing' mesh ");
 635     CPPUNIT_ASSERT((myMesh2->getDescription()).size() != 0);
 636
 637     CPPUNIT_ASSERT(myMesh2->getSpaceDimension() == SpaceDimension);
 638     CPPUNIT_ASSERT(myMesh2->getMeshDimension() == meshDimension);
 639     CPPUNIT_ASSERT(myMesh2->getNumberOfNodes() == NumberOfNodes);
 640
 641     const COORDINATE* coord = myMesh2->getCoordinateptr();
 642     try
 643       {
 644         CPPUNIT_ASSERT(myMesh2->getCoordinatesSystem() != "catresian");
 645       }
 646     catch (const std::exception &e)
 647       {
 648         CPPUNIT_FAIL(e.what());
 649       }
 650     catch (...)
 651       {
 652         CPPUNIT_FAIL("Unknown exception");
 653       }
 654
 655     const string * units;
 656     try
 657       {
 658         units = myMesh2->getCoordinatesUnits();
 659         for (int axe = 0; axe < SpaceDimension; axe++)
 660           {
 661             string verif = coord->getCoordinateUnit(axe+1);
 662             CPPUNIT_ASSERT(verif == units[axe]);
 663           }
 664       }
 665     catch (const std::exception &e)
 666       {
 667         CPPUNIT_FAIL(e.what());
 668       }
 669     catch (...)
 670       {
 671         CPPUNIT_FAIL("Unknown exception");
 672       }
 673
 674     const string * noms;
 675     try
 676       {
 677         noms = myMesh2->getCoordinatesNames();
 678         for (int axe = 0; axe < SpaceDimension; axe++)
 679           {
 680             string verif = coord->getCoordinateName(axe+1);
 681             CPPUNIT_ASSERT(verif == noms[axe]);
 682           }
 683       }
 684     catch (const std::exception &e)
 685       {
 686         CPPUNIT_FAIL(e.what());
 687       }
 688     catch (...)
 689       {
 690         CPPUNIT_FAIL("Unknown exception");
 691       }
 692
 693     try
 694       {
 695         const double * coor2 = myMesh2->getCoordinates(MED_FULL_INTERLACE);
 696
 697         for (int axe = 0; axe < SpaceDimension; axe++)
 698           {
 699             try
 700               {
 701                 for (int num = 0; num < NumberOfNodes; num++)
 702                   {
 703                     try
 704                       {
 705                         const double d = myMesh2->getCoordinate(num + 1, axe + 1);
 706                         CPPUNIT_ASSERT(fabs(d - coor2[(num * SpaceDimension)+axe]) < 0.001);
 707                       }
 708                     catch (const std::exception &e)
 709                       {
 710                         CPPUNIT_FAIL(e.what());
 711                       }
 712                     catch (...)
 713                       {
 714                         CPPUNIT_FAIL("Unknown exception");
 715                       }
 716                   }
 717               }
 718             catch (const std::exception &e)
 719               {
 720                 CPPUNIT_FAIL(e.what());
 721               }
 722             catch (...)
 723               {
 724                 CPPUNIT_FAIL("Unknown exception");
 725               }
 726           }
 727       }
 728     catch (const std::exception &e)
 729       {
 730         CPPUNIT_FAIL(e.what());
 731       }
 732     catch (...)
 733       {
 734         CPPUNIT_FAIL("Unknown exception");
 735       }
 736
 737     const CONNECTIVITY * myConnectivity = myMesh2->getConnectivityptr();
 738
 739     // MED_EN::MED_CELL
 740     MED_EN::medEntityMesh entity = myConnectivity->getEntity();
 741     CPPUNIT_ASSERT_EQUAL(MED_CELL, entity);
 742
 743     int typesNb;
 744     CPPUNIT_ASSERT_NO_THROW(typesNb= myConnectivity->getNumberOfTypes(entity));
 745     CPPUNIT_ASSERT_EQUAL(NumberOfTypes, typesNb);
 746
 747     const MED_EN::medGeometryElement * Types1;
 748     CPPUNIT_ASSERT_NO_THROW(Types1 = myMesh2->getTypes(entity));
 749
 750     medConnectivity myMedConnect;
 751     bool existConnect = false;
 752     if (myMesh2->existConnectivity(MED_NODAL, entity))
 753       {
 754         existConnect = true;
 755         myMedConnect = MED_NODAL;
 756       }
 757     else if(myMesh2->existConnectivity(MED_DESCENDING, entity))
 758       {
 759         existConnect = true;
 760         myMedConnect = MED_DESCENDING;
 761       }
 762
 763     for(int t = 0; t < NumberOfTypes; t++ )
 764       {
 765         CPPUNIT_ASSERT_EQUAL(Types1[t], Types[t]);
 766         int NumberOfElements1 = 0;
 767         CPPUNIT_ASSERT_NO_THROW(NumberOfElements1 = myMesh2->getNumberOfElements(entity, Types1[t]));
 768         CPPUNIT_ASSERT_EQUAL(NumberOfElements1, NumberOfElements[t]);
 769         if(existConnect)
 770           {
 771             const int * connectivity;
 772             const int * connectivity_index;
 773             CPPUNIT_ASSERT_NO_THROW(connectivity = myMesh2->getConnectivity
 774                                     (myMedConnect, entity, Types1[t]));
 775             connectivity_index = myMesh2->getConnectivityIndex(myMedConnect, entity);
 776             for (int j = 0; j < NumberOfElements1; j++)
 777               {
 778                 cout<<"!!!!!!!!!!!!!!!"<<endl;
 779                 for (int k = connectivity_index[j]; k < connectivity_index[j+1]; k++)
 780                   cout << connectivity[k-1] << " ";
 781                 cout << endl;
 782               }
 783           }
 784       }
 785
 786     const CELLMODEL* myCellModel = myMesh2->getCellsTypes(entity);
 787     string* TypeNames;
 788     CPPUNIT_ASSERT_NO_THROW(TypeNames = myMesh2->getCellTypeNames(entity));
 789
 790     for(int k = 0; k < NumberOfTypes; k++ )
 791       {
 792         CPPUNIT_ASSERT_EQUAL(TypeNames[k], myCellModel[k].getName());
 793       }
 794     delete [] TypeNames;
 795
 796     const int* myGlobalNbIdx;
 797     CPPUNIT_ASSERT_NO_THROW(myGlobalNbIdx = myMesh2->getGlobalNumberingIndex(MED_FACE));
 798     for(int i = 0; i <= NumberOfFacesTypes; i++)
 799       {
 800         if(i == NumberOfFacesTypes)
 801           {
 802             CPPUNIT_ASSERT_EQUAL(myGlobalNbIdx[i],nbTria+nbQua+1);
 803             CPPUNIT_ASSERT_THROW(myMesh2->getElementType(MED_FACE, myGlobalNbIdx[i]), MEDEXCEPTION);
 804             break;
 805           }
 806         cout<<"Global number of first element of each geom type : "<<myGlobalNbIdx[i]<<endl;
 807       }
 808     {
 809       const int * ReverseNodalConnectivity;
 810
 811       // Show Reverse Nodal Connectivity
 812       //#ifndef ENABLE_FAULTS
 813       // (BUG) CONNECTIVITY::_numberOfNodes is not set
 814       ((CONNECTIVITY*)myConnectivity)->setNumberOfNodes(NumberOfNodes);
 815       //#endif
 816       //#ifdef ENABLE_FORCED_FAILURES
 817       //      CPPUNIT_FAIL("ERROR in CONNECTIVITY::calculateReverseNodalConnectivity()"
 818       //                   " because myMesh2->_connectivity->_numberOfNodes is not set");
 819       //#endif
 820
 821       CPPUNIT_ASSERT_NO_THROW(ReverseNodalConnectivity = myMesh2->getReverseConnectivity(MED_NODAL, entity));
 822       CPPUNIT_ASSERT_NO_THROW(myMesh2->getReverseConnectivityLength(MED_NODAL, entity));
 823       const int * ReverseNodalConnectivityIndex = myMesh2->getReverseConnectivityIndex(MED_NODAL, entity);
 824       const int ReverseIdxLength = myMesh2->getReverseConnectivityIndexLength(MED_NODAL, entity);
 825       CPPUNIT_ASSERT(ReverseIdxLength == NumberOfNodes+1);
 826       for (int i = 0; i < NumberOfNodes; i++)
 827         {
 828           cout << "Node "<< i+1 << " : ";
 829           for (int j = ReverseNodalConnectivityIndex[i]; j < ReverseNodalConnectivityIndex[i+1]; j++)
 830             cout << ReverseNodalConnectivity[j-1] << " ";
 831           cout << endl;
 832         }
 833
 834       // Show Descending Connectivity
 835       int NumberOfElements1;
 836       const int * connectivity;
 837       const int * connectivity_index;
 838       myMesh2->calculateConnectivity(MED_DESCENDING, entity);
 839       try
 840         {
 841           NumberOfElements1 = myMesh2->getNumberOfElements(entity, MED_ALL_ELEMENTS);
 842           connectivity = myMesh2->getConnectivity( MED_DESCENDING, entity, MED_ALL_ELEMENTS);
 843           connectivity_index = myMesh2->getConnectivityIndex(MED_DESCENDING, entity);
 844         }
 845       catch (MEDEXCEPTION m)
 846         {
 847           CPPUNIT_FAIL(m.what());
 848         }
 849
 850       for (int j = 0; j < NumberOfElements1; j++)
 851         {
 852           cout << "Element " << j+1 << " : ";
 853           for (int k = connectivity_index[j]; k < connectivity_index[j+1]; k++)
 854             cout << connectivity[k-1] << " ";
 855           cout << endl;
 856         }
 857
 858       // getElementNumber
 859       if (myMesh2->existConnectivity(MED_NODAL, MED_FACE))
 860         {
 861           int myTr[3] =
 862             {
 863               1,5,4
 864             };
 865           CPPUNIT_ASSERT_NO_THROW(myMesh2->getElementNumber(MED_NODAL,MED_FACE,MED_TRIA3,myTr));
 866         }
 867     }
 868
 869     //test family and group
 870     int NumberOfGroups;
 871     CPPUNIT_ASSERT_THROW(myMesh2->getNumberOfGroups(MED_ALL_ENTITIES), MEDEXCEPTION);
 872     CPPUNIT_ASSERT_NO_THROW(NumberOfGroups = myMesh2->getNumberOfGroups(MED_CELL));
 873     CPPUNIT_ASSERT_EQUAL(NumberOfGroups, 2);
 874     vector<GROUP*> groups;
 875     CPPUNIT_ASSERT_NO_THROW(groups = myMesh2->getGroups(MED_CELL));
 876     CPPUNIT_ASSERT(groups.size() != 0);
 877     for(int nb = 1; nb <= NumberOfGroups; nb++ )
 878       {
 879         const GROUP* group;
 880         CPPUNIT_ASSERT_NO_THROW(group = myMesh2->getGroup(MED_CELL, nb));
 881         CPPUNIT_ASSERT_EQUAL(group->getName(), groups[nb-1]->getName());
 882       }
 883
 884     int NumberOfFamilies;
 885     CPPUNIT_ASSERT_NO_THROW(NumberOfFamilies = myMesh2->getNumberOfFamilies(MED_CELL));
 886     CPPUNIT_ASSERT_MESSAGE("Current mesh hasn't Families", NumberOfFamilies == 0);
 887
 888     //create families - it's not possible to create, becase not all entities are defined
 889     CPPUNIT_ASSERT_THROW( myMesh2->createFamilies(),MEDEXCEPTION);
 890
 891     /*CPPUNIT_ASSERT_NO_THROW(NumberOfFamilies = myMesh2->getNumberOfFamilies(MED_CELL));
 892       CPPUNIT_ASSERT( NumberOfFamilies != 0);*/
 893
 894     delete myMesh2;
 895   }
 896
 897   //////////////////////////////////////////////////////////////
 898   // TEST 2: Polygon and Polyhedron(only NODAL connectivity)  //
 899   /////////////////////////////////////////////////////////////
 900
 901   double CoordinatesPoly[57] =
 902     {
 903       2.0, 3.0, 2.0,
 904       3.0, 2.0, 2.0,
 905       4.0, 1.0, 2.0,
 906       2.0, 0.0, 2.0,
 907       0.0, 1.0, 2.0,
 908       1.0, 2.0, 2.0,
 909       2.0, 3.0, 1.0,
 910       3.0, 2.0, 0.0,
 911       4.0, 1.0, 0.0,
 912       2.0, 0.0, 0.0,
 913       0.0, 1.0, 0.0,
 914       1.0, 2.0, 0.0,
 915       5.0, 3.0, 2.0,
 916       7.0, 2.0, 2.0,
 917       6.0, 0.0, 2.0,
 918       6.0, 3.0, 0.0,
 919       7.0, 2.0, 0.0,
 920       6.0, 0.0, -1.0,
 921       5.0, 1.0, -3.0
 922     };
 923
 924   const int REFnodalConnOfFaces[91] =
 925     {
 926       1, 2, 3, 4, 5, 6, -1,// Polyhedron 1
 927       1, 7, 8, 2,       -1,
 928       2, 8, 9, 3,       -1,
 929       4, 3, 9, 10,      -1,
 930       5, 4, 10, 11,     -1,
 931       6, 5, 11, 12,     -1,
 932       1, 6, 12, 7,      -1,
 933       7, 12, 8, 10,     -1,
 934       9, 8, 12, 11,
 935
 936       13, 14, 15, 3, 2, -1,// Polyhedron 2
 937       13, 2, 8, 16,     -1,
 938       14, 13, 16, 17,   -1,
 939       15, 14, 17, 15,   -1,
 940       17, 18, 15,       -1,
 941       18, 9, 3,         -1,
 942       15, 9, 2,         -1,
 943       3, 9, 8,          -1,
 944       8, 9, 17, 16,     -1,
 945       9, 18, 17
 946     };
 947   const int NumberOfFaces = 19;
 948   const int NumberOfPolyhedron = 2;
 949   const int nbOfPolygons = 2;
 950
 951   const int REFpolyIndex[NumberOfPolyhedron+1] =
 952     {
 953       1,47,92
 954     };
 955
 956   double PolygonCoordinates[27] =
 957     {
 958       2.0, 3.0, 12.0,
 959       3.0, 2.0, 12.0,
 960       4.0, 1.0, 12.0,
 961       2.0, 0.0, 12.0,
 962       0.0, 1.0, 12.0,
 963       1.0, 2.0, 12.0,
 964       5.0, 3.0, 12.0,
 965       7.0, 2.0, 12.0,
 966       6.0, 0.0, 12.0
 967     };
 968
 969   const int REFpolygonFaces[11] =
 970     {
 971       1, 2, 3, 4, 5, 6, // Polygon 1
 972       7, 8, 9, 3, 2
 973     }; // Polygon 2
 974
 975   const int REFpolygonIndex[nbOfPolygons+1] =
 976     {
 977       1, 7, 12
 978     };
 979
 980   MESHING *myMeshingPoly=new MESHING;
 981   myMeshingPoly->setName("meshingpoly");
 982
 983   const int NbOfTypes = 1;
 984   medGeometryElement TypesPoly[NbOfTypes] =
 985     {
 986       MED_TETRA4
 987     };
 988   const int NbOfElements[NbOfTypes] =
 989     {
 990       1
 991     };
 992
 993   CPPUNIT_ASSERT_NO_THROW(myMeshingPoly->setNumberOfTypes(NbOfTypes, MED_CELL));
 994
 995   try
 996     {
 997       myMeshingPoly->setCoordinates(SpaceDimension, NumberOfNodes, CoordinatesPoly,
 998                                     "CARTESIAN", MED_FULL_INTERLACE);
 999     }
1000   catch (const std::exception &e)
1001     {
1002       CPPUNIT_FAIL(e.what());
1003     }
1004   catch (...)
1005     {
1006       CPPUNIT_FAIL("Unknown exception");
1007     }
1008
1009   CPPUNIT_ASSERT_NO_THROW(myMeshingPoly->setTypes(TypesPoly, MED_CELL));
1010   CPPUNIT_ASSERT_NO_THROW(myMeshingPoly->setNumberOfElements(NbOfElements, MED_CELL));
1011
1012   string Unit ="cm";
1013   for(int i = 0; i < SpaceDimension; i++ )
1014     {
1015       try
1016         {
1017           myMeshingPoly->setCoordinateName(Names[i],i);
1018         }
1019       catch (const std::exception &e)
1020         {
1021           CPPUNIT_FAIL(e.what());
1022         }
1023       catch (...)
1024         {
1025           CPPUNIT_FAIL("Unknown exception");
1026         }
1027       try
1028         {
1029           myMeshingPoly->setCoordinateUnit(Unit, i);
1030         }
1031       catch (const std::exception &e)
1032         {
1033           CPPUNIT_FAIL(e.what());
1034         }
1035       catch (...)
1036         {
1037           CPPUNIT_FAIL("Unknown exception");
1038         }
1039     }
1040
1041   int ConnectivityTetraPoly[4*1]=
1042     {
1043       17, 9, 18, 19
1044     };
1045
1046   CPPUNIT_ASSERT_NO_THROW(myMeshingPoly->setConnectivity( MED_CELL, MED_TETRA4, ConnectivityTetraPoly ));
1047
1048   CPPUNIT_ASSERT_NO_THROW(myMeshingPoly->setConnectivity(MED_CELL, MED_POLYHEDRA,REFnodalConnOfFaces,REFpolyIndex));
1049
1050   bool PolyConn = false;
1051   CPPUNIT_ASSERT_NO_THROW(PolyConn = myMeshingPoly->existConnectivity(MED_NODAL, MED_CELL));
1052   CPPUNIT_ASSERT(PolyConn);
1053     {
1054     CPPUNIT_ASSERT_EQUAL(NumberOfPolyhedron,
1055                          myMeshingPoly->getNumberOfElements(MED_CELL,MED_POLYHEDRA));
1056     CPPUNIT_ASSERT_NO_THROW( myMeshingPoly->calculateConnectivity (MED_NODAL,MED_FACE));
1057     CPPUNIT_ASSERT_EQUAL(NumberOfFaces-1, myMeshingPoly->getNumberOfElements(MED_FACE,MED_POLYGON)); // -1: one face is shared with tetra
1058     CPPUNIT_ASSERT_EQUAL(91,myMeshingPoly->getConnectivityLength(MED_NODAL,MED_CELL,MED_POLYHEDRA));
1059     const int * PolyConn;
1060     const int * PolyIdx;
1061     CPPUNIT_ASSERT_NO_THROW(PolyConn = myMeshingPoly->getConnectivity(MED_NODAL,MED_CELL,MED_ALL_ELEMENTS));
1062     CPPUNIT_ASSERT_NO_THROW(PolyIdx = myMeshingPoly->getConnectivityIndex(MED_NODAL,MED_CELL));
1063     for(int i = NbOfElements[0], iRef=0; i<NbOfElements[0]+NumberOfPolyhedron; i++)
1064       {
1065         int NodeIdxBegin = PolyIdx[i];
1066         int NodeIdxEnd = PolyIdx[i+1];
1067         for(int k = NodeIdxBegin; k < NodeIdxEnd; k++)
1068           CPPUNIT_ASSERT_EQUAL(REFnodalConnOfFaces[iRef++], PolyConn[k-1]);
1069       }
1070   }
1071
1072   MESHING *myPolygonMeshing=new MESHING;
1073   myPolygonMeshing->setName("PolygonMeshing");
1074
1075   medGeometryElement PolygonTypes[NbOfTypes] =
1076     {
1077       MED_TRIA3
1078     };
1079   const int PolygonNumberOfElements[NbOfTypes] =
1080     {
1081       2
1082     };
1083
1084   CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->setNumberOfTypes(NbOfTypes, MED_CELL));
1085
1086   try
1087     {
1088       myPolygonMeshing->setCoordinates(SpaceDimension, NumberOfNodes, PolygonCoordinates,
1089                                        "CARTESIAN", MED_FULL_INTERLACE);
1090     }
1091   catch (const std::exception &e)
1092     {
1093       CPPUNIT_FAIL(e.what());
1094     }
1095   catch (...)
1096     {
1097       CPPUNIT_FAIL("Unknown exception");
1098     }
1099
1100     CPPUNIT_ASSERT_EQUAL(SpaceDimension, myPolygonMeshing->getSpaceDimension());
1101     CPPUNIT_ASSERT_EQUAL(NumberOfNodes, myPolygonMeshing->getNumberOfNodes());
1102
1103     CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->setTypes(PolygonTypes, MED_CELL));
1104     CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->setNumberOfElements(PolygonNumberOfElements, MED_CELL));
1105     CPPUNIT_ASSERT_EQUAL(2, myPolygonMeshing->getMeshDimension());
1106
1107   try
1108     {
1109       myPolygonMeshing->setCoordinatesNames(Names);
1110     }
1111   catch (const std::exception &e)
1112     {
1113       CPPUNIT_FAIL(e.what());
1114     }
1115   catch (...)
1116     {
1117       CPPUNIT_FAIL("Unknown exception");
1118     }
1119
1120   try
1121     {
1122       myPolygonMeshing->setCoordinatesUnits(Units);
1123     }
1124   catch (const std::exception &e)
1125     {
1126       CPPUNIT_FAIL(e.what());
1127     }
1128   catch (...)
1129     {
1130       CPPUNIT_FAIL("Unknown exception");
1131     }
1132
1133   const int sizeTri = 3*2;
1134   int ConnectivityTri[sizeTri]=
1135     {
1136       1, 7, 2, 3, 9, 4
1137     };
1138
1139   CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->setConnectivity( MED_CELL, MED_TRIA3, ConnectivityTri ));
1140   CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->setConnectivity( MED_CELL, MED_POLYGON, REFpolygonFaces, REFpolygonIndex ));
1141
1142   bool PolygonConn = false;
1143   CPPUNIT_ASSERT_NO_THROW(PolygonConn = myPolygonMeshing->existConnectivity(MED_NODAL, MED_CELL));
1144   if(PolygonConn)
1145     {
1146       int Polytypes;
1147       CPPUNIT_ASSERT_NO_THROW(Polytypes = myPolygonMeshing->getNumberOfTypes(MED_CELL));
1148       CPPUNIT_ASSERT_EQUAL(NbOfTypes,Polytypes);
1149
1150       const MED_EN::medGeometryElement * PolyTypes;
1151       CPPUNIT_ASSERT_NO_THROW(PolyTypes = myPolygonMeshing->getTypes(MED_CELL));
1152       CPPUNIT_ASSERT_EQUAL(PolyTypes[NbOfTypes-1],MED_POLYGON);
1153
1154       for(int t = 0; t < Polytypes; t++)
1155         {
1156           CPPUNIT_ASSERT_NO_THROW( myPolygonMeshing->getNumberOfElements(MED_CELL, PolyTypes[t]));
1157         }
1158
1159       medGeometryElement geomPolyElem;
1160       CPPUNIT_ASSERT_NO_THROW(geomPolyElem = myPolygonMeshing->getElementType(MED_CELL, 1));
1161       CPPUNIT_ASSERT_EQUAL(geomPolyElem, MED_TRIA3);
1162
1163       CPPUNIT_ASSERT_EQUAL(myPolygonMeshing->getNumberOfElements(MED_CELL,MED_POLYGON),nbOfPolygons);
1164       CPPUNIT_ASSERT_NO_THROW(myPolygonMeshing->getConnectivityLength(MED_NODAL,MED_CELL,MED_POLYGON));
1165       myPolygonMeshing->removeReference();
1166       const int * PolygonConn;
1167       CPPUNIT_ASSERT_THROW(PolygonConn = myMeshingPoly->getConnectivity(MED_NODAL,MED_CELL,MED_POLYGON),MEDEXCEPTION);
1168     }
1169
1170   ////////////////////////////////////////////////////////////
1171   // TEST : SUPPORT* sup = myMeshPointe->getSupportOnAll()) //
1172   ////////////////////////////////////////////////////////////
1173
1174   //#ifdef ENABLE_FAULTS
1175   {
1176     MESH * myMeshPointe = new MESH(MED_DRIVER, filename, meshname);
1177     const SUPPORT* sup = myMeshPointe->getSupportOnAll(MED_CELL);
1178     CPPUNIT_ASSERT( sup->isOnAllElements() );
1179     CPPUNIT_ASSERT_EQUAL( myMeshPointe->getNumberOfTypes( sup->getEntity() ),
1180                           sup->getNumberOfTypes());
1181     CPPUNIT_ASSERT( sup->getNumber( MED_ALL_ELEMENTS ));
1182     myMeshPointe->removeReference();
1183   }
1184   //#endif
1185   //#ifdef ENABLE_FORCED_FAILURES
1186   //  CPPUNIT_FAIL("ERROR: can not create SUPPORT on mesh, read from pointe.med");
1187   //#endif
1188
1189   ////////////////////////////////////////////////////////
1190   // TEST 3: test MESH on  MEDMEMTest::createTestMesh()//
1191   ///////////////////////////////////////////////////////
1192
1193   MESH* myMesh3 = MEDMEMTest_createTestMesh();
1194
1195   int MeshDim  = myMesh3->getMeshDimension();
1196   medEntityMesh constituentEntity;
1197   if (MeshDim==3)
1198     {
1199       constituentEntity = MED_CELL;
1200     }
1201   if (MeshDim==2)
1202     {
1203       constituentEntity = MED_FACE;
1204     }
1205   if (MeshDim==1)
1206     {
1207       constituentEntity = MED_EDGE;
1208     }
1209
1210   int SpaceDim = myMesh3->getSpaceDimension();
1211
1212   // Show Reverse Nodal Connectivity
1213   const int* ReverseNodalConnectivity;
1214   const int* ReverseNodalConnectivityIndex;
1215   int ReverseLength;
1216   int ReverseIdxLength;
1217
1218   CONNECTIVITY* myConnectivity3 = (CONNECTIVITY*)myMesh3->getConnectivityptr();
1219   myConnectivity3->setNumberOfNodes(myMesh3->getNumberOfNodes());
1220
1221   CPPUNIT_ASSERT_NO_THROW(ReverseNodalConnectivity= myMesh3->getReverseConnectivity(MED_NODAL, MED_CELL));
1222   CPPUNIT_ASSERT_NO_THROW(ReverseLength = myMesh3->getReverseConnectivityLength(MED_NODAL, MED_CELL));
1223   CPPUNIT_ASSERT_NO_THROW(ReverseNodalConnectivityIndex = myMesh3->getReverseConnectivityIndex(MED_NODAL, MED_CELL));
1224   CPPUNIT_ASSERT_NO_THROW(ReverseIdxLength = myMesh3->getReverseConnectivityIndexLength(MED_NODAL, MED_CELL));
1225   CPPUNIT_ASSERT(ReverseIdxLength == myMesh3->getNumberOfNodes()+1);
1226
1227   for (int i = 0; i < myMesh3->getNumberOfNodes(); i++)
1228     {
1229       cout << "Node "<< i+1 << " : ";
1230       for (int j = ReverseNodalConnectivityIndex[i]; j < ReverseNodalConnectivityIndex[i+1]; j++)
1231         cout << ReverseNodalConnectivity[j-1] << " ";
1232       cout << endl;
1233     }
1234
1235   // Show Descending Connectivity
1236   int NumberOfElements1;
1237   const int * connectivity;
1238   const int * connectivity_index;
1239   myMesh3->calculateConnectivity(MED_DESCENDING, MED_EN::MED_CELL);
1240   try
1241     {
1242       NumberOfElements1 = myMesh3->getNumberOfElements(MED_CELL, MED_ALL_ELEMENTS);
1243       connectivity = myMesh3->getConnectivity( MED_DESCENDING, MED_CELL, MED_ALL_ELEMENTS);
1244       connectivity_index = myMesh3->getConnectivityIndex(MED_DESCENDING, MED_CELL);
1245     }
1246   catch (MEDEXCEPTION m)
1247     {
1248       CPPUNIT_FAIL(m.what());
1249     }
1250
1251   for (int j = 0; j < NumberOfElements1; j++)
1252     {
1253       cout << "Element " << j+1 << " : ";
1254       for (int k = connectivity_index[j]; k < connectivity_index[j+1]; k++)
1255         cout << connectivity[k-1] << " ";
1256       cout << endl;
1257     }
1258
1259   //test 3D mesh
1260   for(int ind = SpaceDim; ind > 1; ind-- )
1261     {
1262       int NumberOfElem = myMesh3->getNumberOfElements (constituentEntity,MED_ALL_ELEMENTS);
1263       if(NumberOfElem < 1) continue;
1264
1265       const SUPPORT * sup = myMesh3->getSupportOnAll( constituentEntity );
1266
1267       if (ind == 2)
1268         {
1269           // test of normal(for 1d or 2d elements)
1270           FIELD<double>* normal;
1271           CPPUNIT_ASSERT_NO_THROW(normal = myMesh3->getNormal(sup));
1272
1273           double normal_square, norm;
1274           double maxnorm=0.;
1275           double minnorm=0.;
1276           double tmp_value;
1277           for (int i = 1; i<=NumberOfElem; i++)
1278             {
1279               normal_square = 0.;
1280               cout << "Normal " << i << " ";
1281               for (int j=1; j<=SpaceDim; j++)
1282                 {
1283                   tmp_value = normal->getValueIJ(i,j);
1284                   normal_square += tmp_value*tmp_value;
1285                   cout << tmp_value << " ";
1286                 }
1287               norm = sqrt(normal_square);
1288               maxnorm = dmax(maxnorm,norm);
1289               minnorm = dmin(minnorm,norm);
1290               cout << ", Norm = " << norm << endl;
1291             }
1292           cout << "Max Norm " << maxnorm << " Min Norm " << minnorm << endl;
1293           delete normal;
1294
1295           // test of area(for 2d elements)
1296           FIELD<double>* area;
1297           CPPUNIT_ASSERT_NO_THROW(area = myMesh3->getArea(sup));
1298
1299           double maxarea,minarea,areatot;
1300           maxarea = 0.;
1301           minarea = 0.;
1302           areatot = 0.0;
1303           for (int i = 1; i<=NumberOfElem;i++)
1304             {
1305               cout << "Area " << i << " " << area->getValueIJ(i,1) << endl;
1306               maxarea = dmax(maxarea,area->getValueIJ(i,1));
1307               minarea = dmin(minarea,area->getValueIJ(i,1));
1308               areatot = areatot + area->getValueIJ(i,1);
1309             }
1310
1311           cout << "Max Area " << maxarea << " Min Area " << minarea << endl;
1312           cout << "Support Area " << areatot << endl;
1313
1314           delete area;
1315         }
1316
1317       // test of barycenter(for 3d and 2d elements)
1318       FIELD<double>* barycenter;
1319       CPPUNIT_ASSERT_NO_THROW(barycenter = myMesh3->getBarycenter(sup));
1320
1321       CPPUNIT_ASSERT_NO_THROW(NumberOfElem = myMesh3->getNumberOfElements(constituentEntity,MED_ALL_ELEMENTS));
1322
1323       for (int i = 1; i<=NumberOfElem;i++)
1324         {
1325           if (ind == 3)
1326             cout << "Barycenter " << i << " " << barycenter->getValueIJ(i,1) << " " << barycenter->getValueIJ(i,2) << " " << barycenter->getValueIJ(i,3) << endl;
1327
1328           if (ind == 2)
1329             cout << "Barycenter " << i << " " << barycenter->getValueIJ(i,1) << " " << barycenter->getValueIJ(i,2) << endl;
1330         }
1331       delete barycenter;
1332
1333       // test of volume(for 3d elements)
1334       if (ind == 3)
1335         {
1336           FIELD<double>* volume;
1337           CPPUNIT_ASSERT_NO_THROW(volume= myMesh3->getVolume(sup));
1338
1339           double maxvol,minvol,voltot;
1340           maxvol = 0.;
1341           minvol = 0.;
1342           voltot = 0.0;
1343           for (int i = 1; i<=NumberOfElem;i++)
1344             {
1345               cout << "Volume " << i << " " << volume->getValueIJ(i,1) << endl;
1346               maxvol = dmax(maxvol,volume->getValueIJ(i,1));
1347               minvol = dmin(minvol,volume->getValueIJ(i,1));
1348               voltot = voltot + volume->getValueIJ(i,1);
1349             }
1350
1351           cout << "Max Volume " << maxvol << " Min Volume " << minvol << endl;
1352           cout << "Support Volume " << voltot << endl;
1353
1354           delete volume;
1355
1356           // test of skin
1357           SUPPORT *skin;
1358           CPPUNIT_ASSERT_NO_THROW(skin = myMesh3->getSkin(sup));
1359
1360           //test mergeSupports and intersectSupports. vactor contains only 1 elements
1361           vector<SUPPORT *> myVectSup;
1362           myVectSup.push_back(skin);
1363
1364           //method return a copy of skin object
1365           SUPPORT *copyMergeSkin;
1366           CPPUNIT_ASSERT_NO_THROW(copyMergeSkin = myMesh3->mergeSupports(myVectSup));
1367           try
1368             {
1369               CPPUNIT_ASSERT(copyMergeSkin->deepCompare(*skin));
1370             }
1371           catch (const std::exception &e)
1372             {
1373               CPPUNIT_FAIL(e.what());
1374             }
1375           catch (...)
1376             {
1377               CPPUNIT_FAIL("Unknown exception");
1378             }
1379
1380           //method return a copy of skin object
1381           SUPPORT *copyIntersectSkin;
1382           CPPUNIT_ASSERT_NO_THROW(copyIntersectSkin = myMesh3->intersectSupports(myVectSup));
1383           try
1384             {
1385               CPPUNIT_ASSERT(copyIntersectSkin->deepCompare(*skin));
1386             }
1387           catch (const std::exception &e)
1388             {
1389               CPPUNIT_FAIL(e.what());
1390             }
1391           catch (...)
1392             {
1393               CPPUNIT_FAIL("Unknown exception");
1394             }
1395
1396           skin->removeReference();
1397           copyMergeSkin->removeReference();
1398           copyIntersectSkin->removeReference();
1399         }
1400
1401       constituentEntity++;
1402     }
1403
1404
1405   // Testing length and normal vectors on 1d elements
1406   {
1407     // coordinates
1408     int NumberOfNodes3 = 4;
1409
1410     string Names3[3] =
1411       {
1412         "X","Y","Z"
1413       };
1414     string Units3[3] =
1415       {
1416         "cm","cm","cm"
1417       };
1418
1419     double Coordinates3[4*2] =
1420       {
1421         0.0,  0.0,  // n1
1422         1.0,  1.0,  // n2
1423         0.0,  1.0,  // n3
1424         1.0,  0.0
1425       }; // n4
1426
1427     const int NumberOfEdgeTypes = 1;
1428     MED_EN::medGeometryElement EdgeTypes[NumberOfEdgeTypes] =
1429       {
1430         MED_SEG2
1431       };
1432     const int NumberOfEdges[NumberOfEdgeTypes] =
1433       {
1434         4
1435       };
1436     int ConnectivityEdge[4*2] =
1437       {
1438         1,2, 2,3, 3,4, 4,1
1439       };
1440
1441     // CREATE THE MESH
1442     MEDMEM::MESHING* myMeshing3 = new MEDMEM::MESHING;
1443     myMeshing3->setName("meshing3");
1444     myMeshing3->setCoordinates(/*SpaceDimension*/2, NumberOfNodes3, Coordinates3,
1445                                "CARTESIAN", MED_EN::MED_FULL_INTERLACE);
1446     myMeshing3->setCoordinatesNames(Names3);
1447     myMeshing3->setCoordinatesUnits(Units3);
1448
1449     // define connectivities
1450     //      cell part
1451     const int NumberOfTypes3 = 1;
1452     medGeometryElement Types3[NumberOfTypes3] =
1453       {
1454         MED_QUAD4
1455       };
1456     const int NumberOfElements3[NumberOfTypes3] =
1457       {
1458         1
1459       };
1460
1461     myMeshing3->setNumberOfTypes(NumberOfTypes3,MED_CELL);
1462     myMeshing3->setTypes(Types3,MED_CELL);
1463     myMeshing3->setNumberOfElements(NumberOfElements3,MED_CELL);
1464
1465     int Connectivityquad[1*4] =
1466       {
1467         1,2,3,4
1468       };
1469
1470     myMeshing3->setConnectivity( MED_CELL, MED_QUAD4, Connectivityquad );
1471
1472     myMeshing3->setNumberOfTypes(NumberOfEdgeTypes, MED_EDGE);
1473     myMeshing3->setTypes(EdgeTypes, MED_EDGE);
1474     myMeshing3->setNumberOfElements(NumberOfEdges, MED_EDGE);
1475
1476     myMeshing3->setConnectivity( MED_EDGE, MED_SEG2, ConnectivityEdge );
1477
1478     //test 2D mesh
1479     int NumberOfElem = myMeshing3->getNumberOfElements (MED_EDGE, MED_ALL_ELEMENTS);
1480
1481     const SUPPORT * sup = myMeshing3->getSupportOnAll( MED_EDGE );
1482
1483     // test of normal(for 1d or 2d elements)
1484     FIELD<double>* normal;
1485     CPPUNIT_ASSERT_NO_THROW(normal = myMeshing3->getNormal(sup));
1486
1487     double normal_square, norm;
1488     double maxnorm=0.;
1489     double minnorm=0.;
1490     double tmp_value;
1491     for (int i = 1; i<=NumberOfElem; i++)
1492       {
1493         normal_square = 0.;
1494         cout << "Normal " << i << " ";
1495         for (int j=1; j<=/*SpaceDimension*/2; j++)
1496           {
1497             tmp_value = normal->getValueIJ(i,j);
1498             normal_square += tmp_value*tmp_value;
1499             cout << tmp_value << " ";
1500           }
1501         norm = sqrt(normal_square);
1502         maxnorm = dmax(maxnorm,norm);
1503         minnorm = dmin(minnorm,norm);
1504         cout << ", Norm = " << norm << endl;
1505       }
1506     cout << "Max Norm " << maxnorm << " Min Norm " << minnorm << endl;
1507
1508     // test of length(for 1d elements)
1509     FIELD<double>* length;
1510     CPPUNIT_ASSERT_NO_THROW(length = myMeshing3->getLength(sup));
1511
1512     double length_value,maxlength,minlength;
1513     maxlength = 0;
1514     minlength = 0;
1515     for (int i = 1; i<=NumberOfElem;i++)
1516       {
1517         length_value = length->getValueIJ(i,1);
1518         cout << "Length " << i << " " << length_value << endl;
1519         maxlength = dmax(maxlength,length_value);
1520         minlength = dmin(minlength,length_value);
1521       }
1522     cout << "Max Length " << maxlength << " Min Length " << minlength << endl;
1523
1524     vector< FIELD<double> *> myVectField1;
1525     myVectField1.push_back(normal);
1526     myVectField1.push_back(length);
1527     CPPUNIT_ASSERT_NO_THROW(myMeshing3->mergeFields(myVectField1));
1528
1529     delete normal;
1530     delete length;
1531
1532     //#ifdef ENABLE_FAULTS
1533     {
1534       // (BUG) Segmentation fault if vector is empty
1535       vector<SUPPORT *> myVectSupEmpty;
1536       CPPUNIT_ASSERT_THROW(myMesh3->mergeSupports(myVectSupEmpty), MEDEXCEPTION);
1537     }
1538     //#endif
1539
1540     // test mergeFields method: Fields have the same value type
1541     //intersectSupports and mergeSupports methods
1542     {
1543       SUPPORT * sup1 = new SUPPORT;
1544       sup1->setMesh( myMeshing3 );
1545       sup1->setEntity( MED_EDGE );
1546       SUPPORT * sup2 = new SUPPORT;
1547       sup2->setMesh( myMeshing3 );
1548       sup2->setEntity( MED_EDGE );
1549       MED_EN::medGeometryElement gtEdges[1] =
1550         {
1551           MED_SEG2
1552         };
1553       int nbEdges1[1] =
1554         {
1555           1
1556         };
1557       int edges1[1] =
1558         {
1559           1
1560         };
1561       int nbEdges2[1] =
1562         {
1563           2
1564         };
1565       int edges2[2] =
1566         {
1567           2,3
1568         };
1569       sup1->setpartial("description 1", 1, 1, gtEdges, nbEdges1, edges1);
1570       sup2->setpartial("description 1", 1, 2, gtEdges, nbEdges2, edges2);
1571
1572       vector<SUPPORT *> myVectSup3;
1573       myVectSup3.push_back(sup1);
1574       myVectSup3.push_back(sup2);
1575       //method return a MergeSup on the union of all SUPPORTs in Supports.
1576       SUPPORT *MergeSup;
1577       CPPUNIT_ASSERT_NO_THROW(MergeSup = myMesh3->mergeSupports(myVectSup3));
1578       cout << *MergeSup << endl;
1579       MergeSup->removeReference();
1580
1581       //method return a intersection of all SUPPORTs in IntersectSup
1582       SUPPORT *IntersectSup;
1583       CPPUNIT_ASSERT_NO_THROW(IntersectSup = myMesh3->intersectSupports(myVectSup3));
1584       if (IntersectSup != NULL) cout<< *IntersectSup <<endl;
1585       IntersectSup->removeReference();
1586
1587       FIELD<double> * length1 = myMeshing3->getLength(sup1);
1588       FIELD<double> * length2 = myMeshing3->getLength(sup2);
1589
1590       vector< FIELD<double> *> myVect12;
1591       myVect12.push_back(length1);
1592       myVect12.push_back(length2);
1593
1594       FIELD<double> * length12;
1595       CPPUNIT_ASSERT_NO_THROW(length12 = myMeshing3->mergeFields(myVect12));
1596       delete length12;
1597
1598       sup1->removeReference();
1599       sup2->removeReference();
1600       delete length1;
1601       delete length2;
1602     }
1603   }
1604
1605   /////////////////////////////////////////////////////////
1606   // TEST 4: test MESH constructed from file pointe.med //
1607   ////////////////////////////////////////////////////////
1608   MESH * myMesh4 = new MESH();
1609   myMesh4->setName(meshname);
1610   MED_MESH_RDONLY_DRIVER myMeshDriver (filename, myMesh4);
1611   myMeshDriver.setMeshName(meshname);
1612
1613   //Mesh has no driver->segmentation violation
1614   //CPPUNIT_ASSERT_THROW(myMesh4->read(), MEDEXCEPTION);
1615
1616   //Add an existing MESH driver.
1617   int myDriver4;
1618   CPPUNIT_ASSERT_NO_THROW(myDriver4 = myMesh4->addDriver(myMeshDriver));
1619
1620   //read all objects in the file
1621   CPPUNIT_ASSERT_NO_THROW(myMesh4->read(myDriver4));
1622
1623   if (myMesh4->getIsAGrid())
1624     {
1625       GRID* myGrid = dynamic_cast<GRID*>(myMesh4);
1626       CPPUNIT_ASSERT(myGrid);
1627     }
1628
1629   //myDriver4->DRONLY->can't write
1630   CPPUNIT_ASSERT_THROW(myMesh4->write(myDriver4), MEDEXCEPTION);
1631
1632   // add new driver
1633   int idMeshV21;
1634   CPPUNIT_ASSERT_NO_THROW(idMeshV21 = myMesh4->addDriver(MED_DRIVER,filenameout21));
1635
1636   //Write all the content of the MESH using driver referenced by the integer handler index.
1637   CPPUNIT_ASSERT_NO_THROW(myMesh4->write(idMeshV21));
1638
1639   // remove driver from mesh
1640   CPPUNIT_ASSERT_NO_THROW(myMesh4->rmDriver(myDriver4));
1641   //#ifdef ENABLE_FORCED_FAILURES
1642   //  CPPUNIT_FAIL("ERROR: driver with index idMedV21 has not been removed");
1643   //#endif
1644   // ensure exception is raised on second attempt to remove driver
1645   //CPPUNIT_ASSERT_THROW(myMesh4->rmDriver(myDriver4),MEDEXCEPTION);
1646
1647   // Create a MESH object using a MESH driver of type MED_DRIVER associated with file fileName.
1648   MESH* myMesh5;
1649   CPPUNIT_ASSERT_NO_THROW(myMesh5 = new MESH(MED_DRIVER, filename, meshname));
1650   if(myMesh5->getIsAGrid())
1651     {
1652       GRID* myGrid = dynamic_cast<GRID*>(myMesh4);
1653       CPPUNIT_ASSERT(myGrid);
1654     }
1655
1656   //ensure two meshes constracted from one file in two different ways are equal
1657   CPPUNIT_ASSERT(myMesh5->deepCompare(*myMesh4));
1658
1659   int myDriver6;
1660   MESH* myMesh6 = new MESH();
1661   try
1662     {
1663       myDriver6 = myMesh6->addDriver(MED_DRIVER, filename, meshname, RDONLY);
1664     }
1665   catch (const std::exception &e)
1666     {
1667       CPPUNIT_FAIL(e.what());
1668     }
1669   catch (...)
1670     {
1671       CPPUNIT_FAIL("Unknown exception");
1672     }
1673
1674   try
1675     {
1676       myMesh6->read(myDriver6);
1677     }
1678   catch (const std::exception &e)
1679     {
1680       CPPUNIT_FAIL(e.what());
1681     }
1682   catch (...)
1683     {
1684       CPPUNIT_FAIL("Unknown exception");
1685     }
1686
1687   //ensure two meshes constracted from one file in two different ways are equal
1688   CPPUNIT_ASSERT(myMesh6->deepCompare(*myMesh4));
1689
1690   //test FAMILY
1691   int NumberOfFamilies4;
1692   CPPUNIT_ASSERT_NO_THROW(NumberOfFamilies4 = myMesh6->getNumberOfFamilies(MED_CELL));
1693   CPPUNIT_ASSERT_MESSAGE("Current mesh hasn't Families", NumberOfFamilies4 != 0);
1694
1695   vector<FAMILY*> families4;
1696   CPPUNIT_ASSERT_NO_THROW(families4 = myMesh6->getFamilies(MED_CELL));
1697   CPPUNIT_ASSERT(families4.size() == NumberOfFamilies4);
1698   for(int nb = 1; nb <= NumberOfFamilies4; nb++ )
1699     {
1700       const FAMILY* family;
1701       CPPUNIT_ASSERT_NO_THROW(family = myMesh6->getFamily(MED_CELL, nb));
1702       CPPUNIT_ASSERT_EQUAL(family->getName(), families4[nb-1]->getName());
1703     }
1704
1705   //get support which reference all elements on the boundary of mesh.
1706   SUPPORT*myBndSup;
1707   CPPUNIT_ASSERT_THROW(myMesh6->getBoundaryElements(MED_CELL), MEDEXCEPTION);
1708   //get only face in 3D.
1709   CPPUNIT_ASSERT_NO_THROW(myBndSup = myMesh6->getBoundaryElements(MED_FACE));
1710
1711   //test buildSupportOnElementsFromElementList and buildSupportOnNodeFromElementList
1712   const int * myConnectivityValue6;
1713   CPPUNIT_ASSERT_NO_THROW(myConnectivityValue6 = myMesh6->getReverseConnectivity(MED_DESCENDING));
1714   const int * myConnectivityIndex6;
1715   CPPUNIT_ASSERT_NO_THROW(myConnectivityIndex6 = myMesh6->getReverseConnectivityIndex(MED_DESCENDING));
1716   int numberOfElem6;
1717   CPPUNIT_ASSERT_NO_THROW(numberOfElem6 = myMesh6->getNumberOfElements(MED_FACE,MED_ALL_ELEMENTS));
1718   list<int> myElementsList6;
1719
1720   for (int i=0; i<numberOfElem6; i++)
1721     if (myConnectivityValue6[myConnectivityIndex6[i]] == 0)
1722       {
1723         myElementsList6.push_back(i+1);
1724       }
1725
1726   SUPPORT * mySupportOnNode;
1727   SUPPORT * mySupportOnElem;
1728   CPPUNIT_ASSERT_NO_THROW(mySupportOnElem = myMesh6->buildSupportOnElementsFromElementList(myElementsList6,MED_FACE));
1729   CPPUNIT_ASSERT(mySupportOnElem->deepCompare(*myBndSup));
1730   CPPUNIT_ASSERT_EQUAL(MED_FACE, mySupportOnElem->getEntity());
1731
1732   list<int>::const_iterator iteronelem = myElementsList6.begin();
1733   for (int i = 1; i <= 3; i++, iteronelem++)
1734     {
1735       CPPUNIT_ASSERT_EQUAL(i, mySupportOnElem->getValIndFromGlobalNumber(*iteronelem));
1736     }
1737
1738   CPPUNIT_ASSERT_NO_THROW(mySupportOnNode = myMesh6->buildSupportOnNodeFromElementList(myElementsList6,MED_FACE));
1739   CPPUNIT_ASSERT(mySupportOnNode->deepCompare( *(myMesh6->getBoundaryElements(MED_NODE))));
1740
1741   //sets mesh fields to initial values
1742   try
1743     {
1744       myMesh6->init();
1745     }
1746   catch (const std::exception &e)
1747     {
1748       CPPUNIT_FAIL(e.what());
1749     }
1750   catch (...)
1751     {
1752       CPPUNIT_FAIL("Unknown exception");
1753     }
1754
1755   //ensure two meshes constracted from one file in two different ways are equal
1756   CPPUNIT_ASSERT(!myMesh6->deepCompare(*myMesh4));
1757
1758   //ensure mesh is empty
1759   CPPUNIT_ASSERT(myMesh6->getSpaceDimension() == MED_INVALID);
1760   CPPUNIT_ASSERT(myMesh6->getNumberOfNodes() == MED_INVALID);
1761   CPPUNIT_ASSERT(myMesh6->getCoordinateptr() == NULL);
1762
1763   delete myMesh4;
1764   delete myMesh5;
1765   delete myMesh6;
1766
1767   MESH* myMesh7 = MEDMEMTest_createTestMesh();
1768   vector< vector<double> > myBndBox;
1769   try
1770     {
1771       myBndBox = myMesh7->getBoundingBox();
1772     }
1773   catch (const std::exception &e)
1774     {
1775       CPPUNIT_FAIL(e.what());
1776     }
1777   catch (...)
1778     {
1779       CPPUNIT_FAIL("Unknown exception");
1780     }
1781
1782   cout<<"Bounding box for createTestMesh()"<<endl;
1783   for(int i = 0; i < myBndBox.size(); i++)
1784     {
1785       for(int j = 0; j < myBndBox[i].size(); j++)
1786         cout<<" "<< myBndBox[i][j]<<" ";
1787       cout<<endl;
1788     }
1789
1790   double CoorPoint[3] =
1791     {
1792       0.0,  0.0, 1.0
1793     }; //n2
1794   int idxElem;
1795   try
1796     {
1797       idxElem = myMesh7->getElementContainingPoint(CoorPoint);
1798     }
1799   catch (const std::exception &e)
1800     {
1801       CPPUNIT_FAIL(e.what());
1802     }
1803   catch (...)
1804     {
1805       CPPUNIT_FAIL("Unknown exception");
1806     }
1807   CPPUNIT_ASSERT(idxElem != -1);
1808
1809   double CoorNoPoint[3] =
1810     {
1811       5.0,  0.0, -5.0
1812     }; //there is no such point
1813   int idxNoElem;
1814   try
1815     {
1816       idxNoElem = myMesh7->getElementContainingPoint(CoorNoPoint);
1817     }
1818   catch (const std::exception &e)
1819     {
1820       CPPUNIT_FAIL(e.what());
1821     }
1822   catch (...)
1823     {
1824       CPPUNIT_FAIL("Unknown exception");
1825     }
1826   CPPUNIT_ASSERT(idxNoElem == -1);
1827 }
1828
1829 int main (int argc, char** argv)
1830 {
1831   MEDMEMTest_testMeshAndMeshing();
1832 }