]> SALOME platform Git repositories - tools/solverlab.git/commitdiff
Salome HOME
Improved writing of meshes and fields
authormichael <michael@localhost.localdomain>
Thu, 13 Jan 2022 16:56:21 +0000 (17:56 +0100)
committermichael <michael@localhost.localdomain>
Thu, 13 Jan 2022 16:56:21 +0000 (17:56 +0100)
CDMATH/mesh/inc/Field.hxx
CDMATH/mesh/inc/Mesh.hxx
CDMATH/mesh/src/Field.cxx
CDMATH/mesh/src/Mesh.cxx

index 489588c16c57b5cb8d2f0d96eab25d694082e7c6..fe79b58f2ece0656a9843085bd675a6b5899cfb7 100755 (executable)
@@ -15,7 +15,6 @@ namespace MEDCoupling
   class MEDFileField1TS;
 }
 
-#include "DoubleTab.hxx"
 #include "Vector.hxx"
 #include "Mesh.hxx"
 
@@ -87,7 +86,7 @@ class Field
         *  */
        Field(const std::string meshfileName, EntityType fieldType, 
                  const std::vector<double> Vconstant,const std::string & fieldName = "",
-                  int meshLevel=0, double time=0.0);
+                  int meshLevel=0, double time=0.0, std::string meshName="");
 
     /**
      * constructor with data
@@ -254,7 +253,7 @@ class Field
     /** 
         * \brief Delete the medcoupling mesh to save memory space
      */
-    void deleteMEDCouplingUMesh();
+    void deleteMEDCouplingMesh();
     
     /** 
      * \brief Returns true iff an unstructured mesh has been loaded
@@ -347,7 +346,7 @@ class Field
 
     void writeVTK ( const std::string fileName, bool fromScratch=true ) const ;
 
-    void writeMED ( const std::string fileName, bool fromScratch=true ) const ;
+    void writeMED ( const std::string fileName, bool fromScratch=true ) ;
 
     void writeCSV ( const std::string fileName ) const ;
 
@@ -367,7 +366,6 @@ class Field
        int _numberOfComponents;
        double _time;
        std::string _fieldName;
-       MEDCoupling::MEDFileField1TS *_ff ;//To save the field when the mesh has been deleted
        
     private:
 
index 99ce471b2ace22790187b078a8591e02e7de982a..23eb46e316d53a5f6d13a6948d7c9a703d527363 100644 (file)
@@ -8,9 +8,12 @@
 #ifndef MESH_HXX_
 #define MESH_HXX_
 
+#include <memory>
+
+#include <MCAuto.hxx>
+#include "NormalizedGeometricTypes"
+
 #include "MEDCouplingUMesh.hxx"
-#include "MEDCouplingIMesh.hxx"
-#include "MEDCouplingFieldDouble.hxx"
 
 /**
  * Mesh class is defined by
 
 namespace MEDCoupling
 {
+class MEDFileMesh;
 class MEDFileUMesh;
 class MEDCouplingMesh;
-class MEDCouplingIMesh;
-class MEDCouplingUMesh;
 class DataArrayIdType;
 }
 namespace ParaMEDMEM
 {
 class DataArrayIdType;
 }
-#include <MCAuto.hxx>
-#include "NormalizedGeometricTypes"
-
 class Node;
 class Cell;
 class Face;
@@ -62,12 +61,18 @@ public: //----------------------------------------------------------------
         */
        Mesh ( void ) ;
 
+       /**
+        * \brief constructor with data from a medcoupling mesh
+        * @param medcoupling mesh 
+        */
+       Mesh( MEDCoupling::MCAuto<const MEDCoupling::MEDCouplingMesh> mesh ) ;
+
        /**
         * \brief constructor with data to load a general unstructured mesh
         * @param filename name of mesh file
         * @param meshLevel : relative mesh dimension : 0->cells, 1->Faces etc
         */
-       Mesh ( const std::string filename, int meshLevel=0 ) ;
+       Mesh ( const std::string filename, const std::string & meshName="" , int meshLevel=0) ;
 
        /**
         * \brief constructor with data for a regular 1D grid 
@@ -110,15 +115,12 @@ public: //----------------------------------------------------------------
         */
        Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz, int split_to_tetrahedra_policy=-1, std::string meshName="MESH3D_Regular_Cuboid_Grid") ;
 
-       Mesh( const MEDCoupling::MEDCouplingIMesh* mesh ) ;
-       Mesh( const MEDCoupling::MEDCouplingUMesh* mesh ) ;
-
        /**
         * \brief constructor with data
         * @param filename : file name of mesh med file
         * @param meshLevel : relative mesh dimension : 0->cells, 1->Faces etc
         */
-       void readMeshMed( const std::string filename, int meshLevel=0 ) ;
+       void readMeshMed( const std::string filename, const std::string & meshName="" , int meshLevel=0) ;
 
        /**
         * \brief constructor by copy
@@ -144,29 +146,11 @@ public: //----------------------------------------------------------------
        int getSpaceDimension( void ) const ;
 
        /**
-        * \brief Mesh dimension
+        * \brief return Mesh dimension
         * @return _meshDim
         */
        int getMeshDimension( void ) const ;
 
-       /**
-        * \brief return The nodes in this mesh
-        * @return _nodes
-        */
-       Node* getNodes ( void ) const ;
-
-       /**
-        * \brief return The cells in this mesh
-        * @return _vertices
-        */
-       Cell* getCells ( void ) const ;
-
-       /**
-        * \brief return the faces in this mesh
-        * @return _vertices
-        */
-       Face* getFaces ( void ) const ;
-
        /**
         * \brief return the number of nodes in this mesh
         * @return _numberOfNodes
@@ -239,12 +223,10 @@ public: //----------------------------------------------------------------
 
        double getZMax( void )  const ;
 
-       std::vector<double> getDXYZ() const ;// for structured meshes
-
        std::vector<mcIdType> getCellGridStructure() const;// for structured meshes
 
        /**
-        * \brief surcharge operator =
+        * \brief overload operator =
         * @param mesh : The Mesh object to be copied
         */
        const Mesh& operator= ( const Mesh& mesh ) ;
@@ -267,13 +249,13 @@ public: //----------------------------------------------------------------
 
        /**
         * \brief return the list of face group names
-        * return _faceGroupNames
+        * @return _faceGroupNames
         */
        std::vector<std::string> getNameOfFaceGroups( void )  const ;
 
        /**
         * \brief return the list of node group names
-        * return _nodeGroupNames
+        * @return _nodeGroupNames
         */
        std::vector<std::string> getNameOfNodeGroups( void )  const ;
 
@@ -289,8 +271,8 @@ public: //----------------------------------------------------------------
         */
        std::vector<MEDCoupling::DataArrayIdType *> getNodeGroups( void )  const ;
 
-    /*
-     * Functions to extract boundary nodes and faces Ids
+    /**
+     * \brief Functions to extract boundary nodes and faces Ids
      */
      /**
       *  \brief return the list of boundary faces Ids
@@ -302,8 +284,8 @@ public: //----------------------------------------------------------------
      * @return _boundaryNodeIds
      */
     std::vector< int > getBoundaryNodeIds() const;
-    /*
-     * Functions to extract group nodes and faces ids
+    /**
+     * \brief Functions to extract group nodes and faces ids
      */
      /** 
       * @return list of face group Ids
@@ -322,7 +304,7 @@ public: //----------------------------------------------------------------
        /**
         * \brief write mesh in the MED format
         */
-       void writeMED ( const std::string fileName ) const ;
+       void writeMED ( const std::string fileName, bool fromScratch = true ) const;
 
        void setGroupAtPlan(double value, int direction, double eps, std::string groupName, bool isBoundaryGroup=true) ;
 
@@ -362,7 +344,7 @@ public: //----------------------------------------------------------------
 
     std::vector< std::string > getElementTypesNames() const ;
        /**
-        * \brief Compute the minimum value over all cells of the ratio cell perimeter/cell vaolume
+        * \brief Compute the minimum value over all cells of the ratio cell perimeter/cell volume
         */
     double minRatioVolSurf() const;
     
@@ -374,7 +356,7 @@ public: //----------------------------------------------------------------
     /** 
      * \brief Delete the medcoupling mesh to save memory space
      */
-    void deleteMEDCouplingUMesh();
+    void deleteMEDCouplingMesh();
     
     /** 
      * \brief Returns true iff an unstructured mesh has been loaded
@@ -384,13 +366,14 @@ public: //----------------------------------------------------------------
 private: //----------------------------------------------------------------
 
        MEDCoupling::MEDCouplingUMesh*  setMesh( void ) ;
-       void setGroups( const MEDCoupling::MEDFileUMesh* medmesh, MEDCoupling::MEDCouplingUMesh*  mu) ;//Read all face and node group
+       void setFaceGroups( const MEDCoupling::MEDFileUMesh* medmesh, MEDCoupling::MEDCouplingUMesh*  mu) ;//Read all face groups
+       void setNodeGroups( const MEDCoupling::MEDFileMesh*  medmesh, MEDCoupling::MEDCouplingUMesh*  mu) ;//Read all node groups
        void addNewFaceGroup( const MEDCoupling::MEDCouplingUMesh *m);//adds one face group in the vectors _faceGroups, _faceGroupNames and _faceGroupIds
        
-       /*
-        * The MEDCoupling mesh
+       /**
+        * \brief The MEDCoupling mesh
         */
-       MEDCoupling::MCAuto<MEDCoupling::MEDCouplingMesh> _mesh;
+       MEDCoupling::MCAuto<MEDCoupling::MEDCouplingMesh> _mesh;// This is either a MEDCouplingUMesh or a MEDCouplingStructuredMesh
 
        bool _meshNotDeleted;
        
@@ -406,100 +389,106 @@ private: //----------------------------------------------------------------
         */
        int _meshDim ;
     
-    /*
-     * Structured mesh parameters
+    /**
+     * \brief Signal a structured mesh
+     */
+       bool _isStructured;
+    /**
+     * \brief Number of cells in each direction (Structured meshes)
      */
-
-    bool _isStructured;
-    
-       double _xMin;
-
-       double _xMax;
-
-       double _yMin;
-
-       double _yMax;
-
-       double _zMin;
-
-       double _zMax;
-
        std::vector<mcIdType> _nxyz;
 
-       std::vector<double> _dxyz;
-       /*
-        * The nodes in this mesh.
+       /**
+        * \brief The nodes in this mesh.
         */
-       Node *_nodes;
+       std::shared_ptr<Node> _nodes;
 
-       /*
-        * The number of nodes in this mesh.
+       /**
+        * \brief The number of nodes in this mesh.
         */
        int _numberOfNodes;
 
-       /*
-        * The faces in this mesh.
+       /**
+        * \brief The faces in this mesh.
         */
-       Face *_faces;
+       std::shared_ptr<Face> _faces;
 
-       /*
-        * The numbers of faces in this mesh.
+       /**
+        * \brief The numbers of faces in this mesh.
         */
        int _numberOfFaces;
 
-       /*
-        * The cells in this mesh.
+       /**
+        * \brief The cells in this mesh.
         */
-       Cell *_cells;
+       std::shared_ptr<Cell> _cells;
 
-       /*
-        * The number of cells in this mesh.
+       /**
+        * \brief The number of cells in this mesh.
         */
        int _numberOfCells;
 
-       /*
-        * The number of edges in this mesh.
+       /**
+        * \brief return The nodes in this mesh
+        * @return _nodes
         */
-       int _numberOfEdges;//Useful to deduce the number of non zero coefficients in the finite element matrix 
+       std::shared_ptr<Node> getNodes ( void ) const ;
 
-       /*
-        * The names of face groups.
+       /**
+        * \brief return The cells in this mesh
+        * @return _vertices
+        */
+       std::shared_ptr<Cell> getCells ( void ) const ;
+
+       /**
+        * \brief return the faces in this mesh
+        * @return _vertices
+        */
+       std::shared_ptr<Face> getFaces ( void ) const ;
+
+       /**
+        * \brief The number of edges in this mesh.
+        */
+       int _numberOfEdges;//Useful to deduce the number of non zero coefficients in a finite element matrix 
+
+       /**
+        * \brief The names of face groups.
         */
        std::vector<std::string> _faceGroupNames;
 
-       /*
-        * The names of node groups.
+       /**
+        * \brief The names of node groups.
         */
        std::vector<std::string> _nodeGroupNames;
 
-       /*
-        * The list of face groups.
+       /**
+        * \brief The list of face groups.
         */
        std::vector<MEDCoupling::MEDCouplingUMesh *> _faceGroups;
-       /*
-        * The list of node groups.
+       /**
+        * \brief The list of node groups.
         */
        std::vector<MEDCoupling::DataArrayIdType *> _nodeGroups;
        
-       /*
-        * The list of face id in each face groups.
+       /**
+        * \brief The list of face id in each face groups.
         */
        std::vector< std::vector<int> > _faceGroupsIds;
        
-       /*
-        * The list of node id in each node groups.
+       /**
+        * \brief The list of node id in each node groups.
         */
        std::vector< std::vector<int> > _nodeGroupsIds;
        
-       /*
-        * Elements types (SEG2, TRI3, QUAD4, HEXA6 ...)
+       /**
+        * \brief Elements types (SEG2, TRI3, QUAD4, HEXA6 ...)
         */
        std::vector< INTERP_KERNEL::NormalizedCellType > _eltsTypes;//List of cell types contained in the mesh
        std::vector< std::string > _eltsTypesNames;//List of cell types contained in the mesh
     std::vector< INTERP_KERNEL::NormalizedCellType > getElementTypes() const;    
     
-    /*
-     * Tools to manage periodic boundary conditions in square/cube geometries
+    /**
+     * \brief Tools to manage periodic boundary conditions in square/cube geometries
      */
      bool _indexFacePeriodicSet;
      std::map<int,int> _indexFacePeriodicMap;
index 3252115340379bb654ac2d38b01a44df313cdd5f..2997451a0cb087136bd9513ca7e9273dbecaba43 100755 (executable)
@@ -30,7 +30,6 @@ Field::Field( EntityType typeField )
 //----------------------------------------------------------------------
 {
        _field=NULL;
-       _ff=NULL;
        _typeField=typeField;
        _numberOfComponents=0;
 }
@@ -47,7 +46,6 @@ Field::~Field( void )
 Field::Field(const std::string fieldName, EntityType type, const Mesh& mesh, int numberOfComponents, double time)
 {
        _field = NULL;
-       _ff=NULL;
        _mesh=Mesh(mesh);
        _typeField=type;
        _numberOfComponents=numberOfComponents;
@@ -68,7 +66,7 @@ void Field::buildFieldMemoryStructure()
        }else if(_typeField==NODES)
        {
                _field=MEDCouplingFieldDouble::New(ON_NODES);
-               array->alloc(_mesh.getNumberOfNodes(),_numberOfComponents);
+               array->alloc(mu->getNumberOfNodes(),_numberOfComponents);
                _field->setMesh(mu);
        }else if(_typeField==FACES)
        {
@@ -99,8 +97,6 @@ Field::Field( const std::string filename, EntityType type,
                int iteration, int order, int meshLevel)
 {
        _field = NULL;
-       _ff=NULL;
-       _mesh=Mesh(filename + ".med", meshLevel);
        _typeField=type;
        _fieldName=fieldName;
 
@@ -108,11 +104,10 @@ Field::Field( const std::string filename, EntityType type,
 }
 
 Field::Field(const std::string meshFileName, EntityType type, const std::vector<double> Vconstant, 
-               const std::string & fieldName, int meshLevel, double time )
+               const std::string & fieldName, int meshLevel, double time, std::string meshName )
 {
        _field = NULL;
-       _ff=NULL;
-       _mesh=Mesh(meshFileName + ".med", meshLevel);
+       _mesh=Mesh(meshFileName + ".med", meshName, meshLevel);
        _typeField=type;
        _numberOfComponents=Vconstant.size();
        _time=time;
@@ -130,7 +125,6 @@ Field::Field(const std::string meshFileName, EntityType type, const std::vector<
 Field::Field(const Mesh& M, EntityType type, const Vector Vconstant, const std::string & fieldName, double time)
 {
        _field = NULL;
-       _ff=NULL;
        _mesh=Mesh(M);
        _typeField=type;
        _numberOfComponents=Vconstant.size();
@@ -149,7 +143,6 @@ Field::Field(const Mesh& M, EntityType type, const Vector Vconstant, const std::
 Field::Field(const Mesh& M, EntityType type, const vector<double> Vconstant, const std::string & fieldName, double time) 
 {
        _field = NULL;
-       _ff=NULL;
        _mesh=Mesh(M);
        _typeField=type;
        _numberOfComponents=Vconstant.size();
@@ -172,7 +165,6 @@ Field::Field( int nDim, const vector<double> Vconstant, EntityType type,
                const std::string & fieldName, double time,double epsilon)
 {
        _field = NULL;
-       _ff=NULL;
        _typeField=type;
        _numberOfComponents=Vconstant.size();
        _time=time;
@@ -219,7 +211,6 @@ Field::Field(const Mesh M, const Vector VV_Left, const Vector VV_Right, double d
                throw CdmathException( "Field::Field: Vectors VV_Left and VV_Right have different sizes");
 
        _field = NULL;
-       _ff=NULL;
        _mesh=Mesh(M);
        _typeField=type;
        _numberOfComponents=VV_Left.getNumberOfRows();
@@ -297,7 +288,6 @@ Field::Field(const Mesh M, const Vector Vin, const Vector Vout, double radius,
        }
 
        _field = NULL;
-       _ff=NULL;
        _mesh=Mesh(M);
        _typeField=type;
        _numberOfComponents=Vout.size();
@@ -350,7 +340,6 @@ Field::readFieldMed( const std::string & fileNameRadical,
        size_t iField = 0;
        std::string attributedFieldName;
        _field = NULL;
-       _ff=NULL;
 
        // Get the name of the right field that we will attribute to the Field.
        if (fieldName == "") {
@@ -430,7 +419,8 @@ Field::readFieldMed( const std::string & fileNameRadical,
        _numberOfComponents = _field->getNumberOfComponents() ;
        _time = _field->getTime(iteration, order);
 
-       cout<<"Found field " << fieldName << " in file " << completeFileName <<endl;
+       cout<<"Found field " << fieldName << " in file " << completeFileName << " on mesh named "<< _field->getMesh()->getName()<< endl;
+       _mesh=Mesh( completeFileName, _field->getMesh()->getName());
 }
 
 
@@ -992,9 +982,6 @@ void
 Field::writeVTK (std::string fileName, bool fromScratch) const
 //----------------------------------------------------------------------
 {
-       if( !_mesh.isStructured() && !_mesh.meshNotDeleted() )
-               throw CdmathException("Field::writeVTK : Cannot save field in VTK format : unstructured mesh with no MEDCouplingUMesh loaded. Use med format.");
-
        string fname=fileName+".pvd";
        int iter,order;
        double time=_field->getTime(iter,order);
@@ -1128,7 +1115,7 @@ Field::writeCSV ( const std::string fileName ) const
 
 //----------------------------------------------------------------------
 void
-Field::writeMED ( const std::string fileName, bool fromScratch) const
+Field::writeMED ( const std::string fileName, bool fromScratch)
 //----------------------------------------------------------------------
 {
        string fname=fileName+".med";
@@ -1138,20 +1125,16 @@ Field::writeMED ( const std::string fileName, bool fromScratch) const
                        MEDCoupling::WriteField(fname.c_str(),_field,fromScratch);
                else
                        MEDCoupling::WriteFieldUsingAlreadyWrittenMesh(fname.c_str(),_field);
-       else//The mesh has ben deleted, use _ff instead of _field to save the values
+       else//The mesh has ben deleted, use a MEDFileField1TS instead of _field to save the values
        {
-               //MEDFileUMesh * meshMEDFile = MEDFileUMesh::New();
-               //meshMEDFile->setMeshAtLevel(0,_field->getMesh()->buildUnstructured());
-               //meshMEDFile->write(fname.c_str(), fromScratch);
-               //MEDCoupling::WriteUMesh(fname.c_str(),_field->getMesh()->buildUnstructured(),fromScratch);
-               //MEDCoupling::WriteMesh(fname.c_str(),_field->getMesh(),fromScratch);
-           //MEDCoupling::MEDCouplingUMesh* fmesh = dynamic_cast<MEDCoupling::MEDCouplingUMesh*> (_field->getMesh()->deepCopy());
-               //cout<<" checkConsecutiveCellTypes : "<< fmesh->checkConsecutiveCellTypes() <<endl;
-               //cout<<" advancedRepr() : "<< fmesh->advancedRepr() <<endl;
-               //cout<<" checkConsecutiveCellTypes : "<< _field->getMesh()->buildUnstructured()->checkConsecutiveCellTypes()<<endl;
-               MEDFileField1TS *ff=MEDFileField1TS::New();
-               _ff->setFieldNoProfileSBT(  _field );
-               _ff->write(fname.c_str(), fromScratch);
+               if ( not fromScratch)
+               {
+                       MEDCoupling::MCAuto<MEDCoupling::MEDFileField1TS> ff=MEDFileField1TS::New();//To save the field when the mesh has been deleted
+                       ff->setFieldNoProfileSBT(  _field );
+                       ff->write(fname.c_str(), fromScratch);
+               }
+               else
+                       throw CdmathException("Field::writeMED Error !!! The mesh has been deleted, cannot write field from scratch");
        }
 }
 
@@ -1224,7 +1207,7 @@ std::ostream& operator<<(std::ostream& out, const Field& field )
        return out;
 }
 
-void Field::deleteMEDCouplingUMesh()
+void Field::deleteMEDCouplingMesh()
 { 
-       return _mesh.deleteMEDCouplingUMesh();
+       return _mesh.deleteMEDCouplingMesh();
 }
index 44ca0405d8b8c3bd2f524ca884452a2915bc3040..8976721a5191a446e0c97c239b35ab0757246c66 100644 (file)
@@ -5,22 +5,26 @@
  *      Authors: CDMATH
  */
 
+#include <iostream>
+#include <stdio.h>
+#include <stdlib.h>
+#include <iterator>
+#include <algorithm> 
+#include <cassert> 
+
 #include "Mesh.hxx"
 #include "Node.hxx"
 #include "Cell.hxx"
 #include "Face.hxx"
 
-#include "CdmathException.hxx"
+#include "MEDCouplingIMesh.hxx"
+#include "MEDCouplingUMesh.hxx"
+#include "MEDCouplingFieldDouble.hxx"
+
 #include "MEDFileMesh.hxx"
 #include "MEDLoader.hxx"
 
-#include <iostream>
-#include <stdio.h>
-#include <stdlib.h>
-#include <iterator>
-#include <algorithm> 
-#include <string.h>
-#include <assert.h>
+#include "CdmathException.hxx"
 
 using namespace MEDCoupling;
 using namespace std;
@@ -41,14 +45,7 @@ Mesh::Mesh( void )
        _numberOfCells = 0;
        _numberOfEdges = 0;
     _isStructured=false;
-       _xMin=0.;
-       _xMax=0.;
-       _yMin=0.;
-       _yMax=0.;
-       _zMin=0.;
-       _zMax=0.;
     _nxyz.resize(0);
-    _dxyz.resize(0.);
        _boundaryMesh=NULL;
     _boundaryFaceIds.resize(0);
     _boundaryNodeIds.resize(0);
@@ -67,10 +64,6 @@ Mesh::Mesh( void )
 Mesh::~Mesh( void )
 //----------------------------------------------------------------------
 {
-       delete [] _cells;
-       delete [] _nodes;
-       delete [] _faces;
-       
        //for(int i=0; i< _faceGroups.size(); i++)
        //      _faceGroups[i]->decrRef();
        //      _nodeGroups[i]->decrRef();
@@ -84,97 +77,34 @@ Mesh::getName( void ) const
     return _name;
 }
 
-Mesh::Mesh( const MEDCoupling::MEDCouplingIMesh* mesh )
-{
-       _spaceDim=mesh->getSpaceDimension();
-       _meshDim=mesh->getMeshDimension();
-       _dxyz=mesh->getDXYZ();
-       _nxyz=mesh->getCellGridStructure();
-       double* Box0=new double[2*_spaceDim];
-       mesh->getBoundingBox(Box0);
-    _name=mesh->getName();
-    _epsilon=1e-6;
-    _indexFacePeriodicSet=false;
-    
-       _xMin=Box0[0];
-       _xMax=Box0[1];
-       if (_spaceDim>=2)
-       {
-               _yMin=Box0[2];
-               _yMax=Box0[3];
-       }
-       if (_spaceDim>=3)
-       {
-               _zMin=Box0[4];
-               _zMax=Box0[5];
-       }
-
-       double *originPtr = new double[_spaceDim];
-       double *dxyzPtr = new double[_spaceDim];
-       mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
-
-       for(int i=0;i<_spaceDim;i++)
-       {
-               originPtr[i]=Box0[2*i];
-               nodeStrctPtr[i]=_nxyz[i]+1;
-               dxyzPtr[i]=_dxyz[i];
-       }
-       _mesh=MEDCouplingIMesh::New(_name,
-                       _spaceDim,
-                       nodeStrctPtr,
-                       nodeStrctPtr+_spaceDim,
-                       originPtr,
-                       originPtr+_spaceDim,
-                       dxyzPtr,
-                       dxyzPtr+_spaceDim);
-    _meshNotDeleted=true;
-    _isStructured=true;
-
-       delete [] originPtr;
-       delete [] dxyzPtr;
-       delete [] nodeStrctPtr;
-       delete [] Box0 ;
-
-       setMesh();
-}
-
-Mesh::Mesh( const MEDCoupling::MEDCouplingUMesh* mesh )
+Mesh::Mesh( MEDCoupling::MCAuto<const MEDCoupling::MEDCouplingMesh> mesh )
 {
        _spaceDim=mesh->getSpaceDimension();
        _meshDim=mesh->getMeshDimension();
-       double* Box0=new double[2*_spaceDim];
-       mesh->getBoundingBox(Box0);
     _name=mesh->getName();
     _epsilon=1e-6;
     _indexFacePeriodicSet=false;
+       _meshNotDeleted=true;
     
-       _xMin=Box0[0];
-       _xMax=Box0[1];
-       if (_spaceDim>=2)
-       {
-               _yMin=Box0[2];
-               _yMax=Box0[3];
-       }
-       if (_spaceDim>=3)
-       {
-               _zMin=Box0[4];
-               _zMax=Box0[5];
-       }
+       _mesh= mesh->clone(false);//Clone because you will need to buildUnstructured. No deep copy : it is assumed node coordinates and cell connectivity will not change
 
-       _mesh=mesh->deepCopy();
-       _mesh=mesh->buildUnstructured();
-    _meshNotDeleted=true;
-    _isStructured=false;
-       delete [] Box0 ;
+    MEDCoupling::MEDCouplingStructuredMesh* structuredMesh = dynamic_cast<MEDCoupling::MEDCouplingStructuredMesh*> (_mesh.retn());
+    if(structuredMesh)
+    {
+        _isStructured=true;
+        _nxyz=structuredMesh->getCellGridStructure();
+    }
+    else
+        _isStructured=false;
 
        setMesh();
 }
 
 //----------------------------------------------------------------------
-Mesh::Mesh( const std::string filename, int meshLevel )
+Mesh::Mesh( const std::string filename, const std::string & meshName, int meshLevel)
 //----------------------------------------------------------------------
 {
-       readMeshMed(filename, meshLevel);
+       readMeshMed(filename, meshName, meshLevel);
 }
 
 //----------------------------------------------------------------------
@@ -185,17 +115,9 @@ Mesh::Mesh( const Mesh& mesh )
        _meshDim = mesh.getMeshDimension() ;
     _name=mesh.getName();
     _epsilon=mesh.getComparisonEpsilon();
-    _xMax=mesh.getXMax();
-    _yMin=mesh.getYMin();
-    _yMax=mesh.getYMax();
-    _zMin=mesh.getZMin();
-    _zMax=mesh.getZMax();
     _isStructured=mesh.isStructured();
     if(_isStructured)
-    {
         _nxyz = mesh.getCellGridStructure() ;
-        _dxyz=mesh.getDXYZ();
-    }
        _numberOfNodes = mesh.getNumberOfNodes();
        _numberOfFaces = mesh.getNumberOfFaces();
        _numberOfCells = mesh.getNumberOfCells();
@@ -206,23 +128,10 @@ Mesh::Mesh( const Mesh& mesh )
        _nodeGroupNames = mesh.getNameOfNodeGroups() ;
        _nodeGroups = mesh.getNodeGroups() ;
 
-       _nodes   = new Node[_numberOfNodes] ;
-       _faces   = new Face[_numberOfFaces] ;
-       _cells   = new Cell[_numberOfCells] ;
+       _nodes   = mesh.getNodes() ;
+       _faces   = mesh.getFaces() ;
+       _cells   = mesh.getCells() ;
     
-    //memcpy(_nodes,mesh.getNodes(),sizeof(*mesh.getNodes())) ;
-    //memcpy(_cells,mesh.getCells(),sizeof(*mesh.getCells())) ;
-    //memcpy(_faces,mesh.getFaces(),sizeof(*mesh.getFaces())) ;
-
-       for (int i=0;i<_numberOfNodes;i++)
-               _nodes[i]=mesh.getNode(i);
-
-       for (int i=0;i<_numberOfFaces;i++)
-               _faces[i]=mesh.getFace(i);
-
-       for (int i=0;i<_numberOfCells;i++)
-               _cells[i]=mesh.getCell(i);
-
     _indexFacePeriodicSet= mesh.isIndexFacePeriodicSet();
     if(_indexFacePeriodicSet)
         _indexFacePeriodicMap=mesh.getIndexFacePeriodic();
@@ -233,17 +142,23 @@ Mesh::Mesh( const Mesh& mesh )
     _eltsTypes=mesh.getElementTypes();
     _eltsTypesNames=mesh.getElementTypesNames();
     
-       MCAuto<MEDCouplingMesh> m1=mesh.getMEDCouplingMesh()->deepCopy();
+       MCAuto<MEDCouplingMesh> m1=mesh.getMEDCouplingMesh()->clone(false);//Clone because you will need to buildUnstructured. No deep copy : it is assumed node coordinates and cell connectivity will not change
+
        _mesh=m1;
     _meshNotDeleted=mesh.meshNotDeleted();
 }
 
 //----------------------------------------------------------------------
 void
-Mesh::readMeshMed( const std::string filename, const int meshLevel)
+Mesh::readMeshMed( const std::string filename, const std::string & meshName, int meshLevel)
 //----------------------------------------------------------------------
 {
-       MEDFileUMesh *m=MEDFileUMesh::New(filename.c_str());//reads the first mesh encountered in the file, otherwise call New (const char *fileName, const char *mName, int dt=-1, int it=-1)
+       MEDFileMesh *m;
+       if( meshName == "" )
+               m=MEDFileMesh::New(filename.c_str());//reads the first mesh encountered in the file, otherwise call New (const char *fileName, const char *mName, int dt=-1, int it=-1)
+       else
+               m=MEDFileMesh::New(filename.c_str(), meshName.c_str());//seeks the mesh named meshName in the file
+
        _mesh=m->getMeshAtLevel(meshLevel);
     _mesh->checkConsistencyLight();
        _mesh->setName(_mesh->getName());
@@ -252,37 +167,22 @@ Mesh::readMeshMed( const std::string filename, const int meshLevel)
     _name=_mesh->getName();
     _epsilon=1e-6;
     _indexFacePeriodicSet=false;
-    MEDCoupling::MEDCouplingIMesh* structuredMesh = dynamic_cast<MEDCoupling::MEDCouplingIMesh*> (_mesh.retn());
+       _meshNotDeleted=true;
+
+    MEDCoupling::MEDCouplingStructuredMesh* structuredMesh = dynamic_cast<MEDCoupling::MEDCouplingStructuredMesh*> (_mesh.retn());
     if(structuredMesh)
     {
         _isStructured=true;
-               _meshNotDeleted=false;
-        _dxyz=structuredMesh->getDXYZ();
         _nxyz=structuredMesh->getCellGridStructure();
-        double* Box0=new double[2*_spaceDim];
-        structuredMesh->getBoundingBox(Box0);
-    
-        _xMin=Box0[0];
-        _xMax=Box0[1];
-        if (_spaceDim>=2)
-        {
-            _yMin=Box0[2];
-            _yMax=Box0[3];
-        }
-        if (_spaceDim>=3)
-        {
-            _zMin=Box0[4];
-            _zMax=Box0[5];
-        }
     }
     else
-    {
         _isStructured=false;
-               _meshNotDeleted=true;
-    }
     
        MEDCouplingUMesh*  mu = setMesh();
-       setGroups(m, mu);
+       setNodeGroups(m, mu);//Works for both cartesan and unstructured meshes
+       MEDFileUMesh *umedfile=dynamic_cast< MEDFileUMesh * > (m);
+       if(umedfile)
+               setFaceGroups(umedfile, mu);//Works only for unstructured meshes
 
        cout<<endl<< "Loaded file "<< filename<<endl;
     cout<<"Mesh name = "<<m->getName()<<", mesh dim = "<< _meshDim<< ", space dim = "<< _spaceDim<< ", nb cells= "<<getNumberOfCells()<< ", nb nodes= "<<getNumberOfNodes()<<endl;
@@ -290,523 +190,231 @@ Mesh::readMeshMed( const std::string filename, const int meshLevel)
        m->decrRef();
 }
 
-void
-Mesh::setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup)
+//----------------------------------------------------------------------
+Mesh::Mesh( std::vector<double> points, std::string meshName )
+//----------------------------------------------------------------------
 {
-       std::vector< int > faceIds(0);
-       double FX, FY, FZ;
-       Face Fi;
-       
-       /* Construction of the face group */
-       if(isBoundaryGroup)        
-        for(int i=0; i<_boundaryFaceIds.size(); i++)
-        {
-                       Fi=_faces[_boundaryFaceIds[i]];
-                       FX=Fi.x();
-                       FY=Fi.y();
-                       FZ=Fi.z();
-                       if (abs(FX-x)<eps && abs(FY-y)<eps && abs(FZ-z)<eps)
-                       {
-                               faceIds.insert(faceIds.end(),_boundaryFaceIds[i]);
-                               _faces[_boundaryFaceIds[i]].setGroupName(groupName);
-                       }
-        }
-       else
-               for (int iface=0;iface<_numberOfFaces;iface++)
-               {
-                       Fi=_faces[iface];
-                       FX=Fi.x();
-                       FY=Fi.y();
-                       FZ=Fi.z();
-                       if (abs(FX-x)<eps && abs(FY-y)<eps && abs(FZ-z)<eps)
-                       {
-                               faceIds.insert(faceIds.end(),iface);
-                               _faces[iface].setGroupName(groupName);
-                       }
-               }
-
-       if (faceIds.size()>0)
+    int nx=points.size();
+    
+       if(nx<2)
+               throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : nx < 2, vector should contain at least two values");
+    int i=0;
+    while( i<nx-1 && points[i+1]>points[i] )
+        i++;
+       if( i!=nx-1 )
     {
-               std::vector< std::string >::iterator it = std::find(_faceGroupNames.begin(), _faceGroupNames.end(), groupName);
-               if(it == _faceGroupNames.end())//No group named groupName
-               {
-                       _faceGroupNames.insert(_faceGroupNames.end(),groupName);
-                       _faceGroupsIds.insert(  _faceGroupsIds.end(),faceIds);
-                       _faceGroups.insert(    _faceGroups.end(), NULL);//No mesh created. Create one ?
-               }
-               else
-               {
-                       std::vector< int > faceGroupIds = _faceGroupsIds[it-_faceGroupNames.begin()];
-                       faceGroupIds.insert( faceGroupIds.end(), faceIds.begin(), faceIds.end());
-                       /* Detect and erase duplicates face ids */
-                       sort( faceGroupIds.begin(), faceGroupIds.end() );
-                       faceGroupIds.erase( unique( faceGroupIds.begin(), faceGroupIds.end() ), faceGroupIds.end() );
-                       _faceGroupsIds[it-_faceGroupNames.begin()] = faceGroupIds;
-               }
-       }
+        //cout<< points << endl;
+               throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : vector values should be sorted");
+    }
+    
+       _spaceDim = 1 ;
+       _meshDim  = 1 ;
+    _name=meshName;
+    _epsilon=1e-6;
+    _indexFacePeriodicSet=false;
+    
+    MEDCouplingUMesh * mesh1d = MEDCouplingUMesh::New(meshName, 1);
+    mesh1d->allocateCells(nx - 1);
+    double * coords = new double[nx];
+    mcIdType nodal_con[2];
+    coords[0]=points[0];
+    for(int i=0; i<nx- 1 ; i++)
+    {
+        nodal_con[0]=i;
+        nodal_con[1]=i+1;
+        mesh1d->insertNextCell(INTERP_KERNEL::NORM_SEG2, 2, nodal_con);
+        coords[i+1]=points[i + 1];
+    }
+    mesh1d->finishInsertingCells();
+
+    DataArrayDouble * coords_arr = DataArrayDouble::New();
+    coords_arr->alloc(nx,1);
+    std::copy(coords,coords+nx,coords_arr->getPointer());
+    mesh1d->setCoords(coords_arr);
+
+    delete [] coords;
+    coords_arr->decrRef();
+
+    _mesh=mesh1d;//To enable writeMED. Because we declared the mesh as unstructured, we decide to build the unstructured data (not mandatory)
+    _meshNotDeleted=true;
+    _isStructured = false;
+
+       setMesh();
 }
 
-void
-Mesh::setGroupAtNodeByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup)
+//----------------------------------------------------------------------
+Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName )
+//----------------------------------------------------------------------
 {
-       std::vector< int > nodeIds(0);
-       double NX, NY, NZ;
-       Node Ni;
-       
-       /* Construction of the node group */
-       if(isBoundaryGroup)        
-        for(int i=0; i<_boundaryNodeIds.size(); i++)
-        {
-                       Ni=_nodes[_boundaryNodeIds[i]];
-                       NX=Ni.x();
-                       NY=Ni.y();
-                       NZ=Ni.z();
-                       if (abs(NX-x)<eps && abs(NY-y)<eps && abs(NZ-z)<eps)
-                       {
-                               nodeIds.insert(nodeIds.end(),_boundaryNodeIds[i]);
-                               _nodes[_boundaryNodeIds[i]].setGroupName(groupName);
-                       }
-        }
-       else
-               for (int inode=0;inode<_numberOfNodes;inode++)
-               {
-                       NX=_nodes[inode].x();
-                       NY=_nodes[inode].y();
-                       NZ=_nodes[inode].z();
-                       if (abs(NX-x)<eps && abs(NY-y)<eps && abs(NZ-z)<eps)
-                       {
-                               nodeIds.insert(nodeIds.end(),inode);
-                               _nodes[inode].setGroupName(groupName);
-                       }
-               }
+       if(nx<=0)
+               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : nx <= 0");
+       if(xmin>=xmax)
+               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : xmin >= xmax");
 
-       if (nodeIds.size()>0)
-    {
-               std::vector< std::string >::iterator it = std::find(_nodeGroupNames.begin(), _nodeGroupNames.end(), groupName);
-               if(it == _nodeGroupNames.end())//No group named groupName
-               {
-                       _nodeGroupNames.insert(_nodeGroupNames.end(),groupName);
-                       _nodeGroupsIds.insert(  _nodeGroupsIds.end(),nodeIds);
-                       _nodeGroups.insert(    _nodeGroups.end(), NULL);//No mesh created. Create one ?
-               }
-               else
-               {
-                       std::vector< int > nodeGroupIds = _nodeGroupsIds[it-_nodeGroupNames.begin()];
-                       nodeGroupIds.insert( nodeGroupIds.end(), nodeIds.begin(), nodeIds.end());
-                       /* Detect and erase duplicates node ids */
-                       sort( nodeGroupIds.begin(), nodeGroupIds.end() );
-                       nodeGroupIds.erase( unique( nodeGroupIds.begin(), nodeGroupIds.end() ), nodeGroupIds.end() );
-                       _nodeGroupsIds[it-_nodeGroupNames.begin()] = nodeGroupIds;
-               }
-    }
+       double dx = (xmax - xmin)/nx ;
+
+       _spaceDim = 1 ;
+       _meshDim  = 1 ;
+    _name=meshName;
+    _epsilon=1e-6;
+    _indexFacePeriodicSet=false;
+
+       _nxyz.resize(_spaceDim);
+       _nxyz[0]=nx;
+
+       double originPtr[_spaceDim];
+       double dxyzPtr[_spaceDim];
+       mcIdType nodeStrctPtr[_spaceDim];
+
+       originPtr[0]=xmin;
+       nodeStrctPtr[0]=nx+1;
+       dxyzPtr[0]=dx;
+
+       _mesh=MEDCouplingIMesh::New(meshName,
+                       _spaceDim,
+                       nodeStrctPtr,
+                       nodeStrctPtr+_spaceDim,
+                       originPtr,
+                       originPtr+_spaceDim,
+                       dxyzPtr,
+                       dxyzPtr+_spaceDim);
+    _meshNotDeleted=true;
+    _isStructured = true;
+
+       setMesh();
 }
 
-void
-Mesh::setGroupAtPlan(double value, int direction, double eps, std::string groupName, bool isBoundaryGroup)
+//----------------------------------------------------------------------
+Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, int split_to_triangles_policy, std::string meshName)
+//----------------------------------------------------------------------
 {
-       std::vector< int > faceIds(0), nodeIds(0);
-       double cord;
-       
-       /* Construction of the face group */    
-       if(isBoundaryGroup)        
-        for(int i=0; i<_boundaryFaceIds.size(); i++)
-        {
-                       cord=_faces[_boundaryFaceIds[i]].getBarryCenter()[direction];
-                       if (abs(cord-value)<eps)
-                       {
-                               faceIds.insert(faceIds.end(),_boundaryFaceIds[i]);
-                               _faces[_boundaryFaceIds[i]].setGroupName(groupName);
-                       }
-        }
-       else
-               for (int iface=0;iface<_numberOfFaces;iface++)
-               {
-                       cord=_faces[iface].getBarryCenter()[direction];
-                       if (abs(cord-value)<eps)
-                       {
-                               faceIds.insert(faceIds.end(),iface);
-                               _faces[iface].setGroupName(groupName);
-                       }
-               }
+       if(nx<=0 || ny<=0)
+               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : nx <= 0 or ny <= 0");
+       if(xmin>=xmax)
+               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : xmin >= xmax");
+       if(ymin>=ymax)
+               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : ymin >= ymax");
 
-       /* Construction of the node group */
-       if(isBoundaryGroup)        
-        for(int i=0; i<_boundaryNodeIds.size(); i++)
-        {
-                       cord=_nodes[_boundaryNodeIds[i]].getPoint()[direction];
-                       if (abs(cord-value)<eps)
-                       {
-                               nodeIds.insert(nodeIds.end(),_boundaryNodeIds[i]);
-                               _nodes[_boundaryNodeIds[i]].setGroupName(groupName);
-                       }
-        }
-       else
-               for (int inode=0;inode<_numberOfNodes;inode++)
-               {
-                       cord=_nodes[inode].getPoint()[direction];
-                       if (abs(cord-value)<eps)
-                       {
-                               nodeIds.insert(nodeIds.end(),inode);
-                               _nodes[inode].setGroupName(groupName);
-                       }
-               }
+       double dx = (xmax - xmin)/nx ;
+       double dy = (ymax - ymin)/ny ;
 
-       if (faceIds.size()>0)
-    {
-               std::vector< std::string >::iterator it = std::find(_faceGroupNames.begin(), _faceGroupNames.end(), groupName);
-               if(it == _faceGroupNames.end())//No group named groupName
-               {
-                       _faceGroupNames.insert(_faceGroupNames.end(),groupName);
-                       _faceGroupsIds.insert(  _faceGroupsIds.end(),faceIds);
-                       _faceGroups.insert(    _faceGroups.end(), NULL);//No mesh created. Create one ?
-               }
-               else
-               {
-                       std::vector< int > faceGroupIds = _faceGroupsIds[it-_faceGroupNames.begin()];
-                       faceGroupIds.insert( faceGroupIds.end(), faceIds.begin(), faceIds.end());
-                       /* Detect and erase duplicates face ids */
-                       sort( faceGroupIds.begin(), faceGroupIds.end() );
-                       faceGroupIds.erase( unique( faceGroupIds.begin(), faceGroupIds.end() ), faceGroupIds.end() );
-                       _faceGroupsIds[it-_faceGroupNames.begin()] = faceGroupIds;
-               }
-       }
-       if (nodeIds.size()>0)
-    {
-               std::vector< std::string >::iterator it = std::find(_nodeGroupNames.begin(), _nodeGroupNames.end(), groupName);
-               if(it == _nodeGroupNames.end())//No group named groupName
-               {
-                       _nodeGroupNames.insert(_nodeGroupNames.end(),groupName);
-                       _nodeGroupsIds.insert(  _nodeGroupsIds.end(),nodeIds);
-                       _nodeGroups.insert(    _nodeGroups.end(), NULL);//No mesh created. Create one ?
-               }
-               else
-               {
-                       std::vector< int > nodeGroupIds = _nodeGroupsIds[it-_nodeGroupNames.begin()];
-                       nodeGroupIds.insert( nodeGroupIds.end(), nodeIds.begin(), nodeIds.end());
-                       /* Detect and erase duplicates node ids */
-                       sort( nodeGroupIds.begin(), nodeGroupIds.end() );
-                       nodeGroupIds.erase( unique( nodeGroupIds.begin(), nodeGroupIds.end() ), nodeGroupIds.end() );
-                       _nodeGroupsIds[it-_nodeGroupNames.begin()] = nodeGroupIds;
-               }
-    }
-}
+       _spaceDim = 2 ;
+       _meshDim  = 2 ;
+    _name=meshName;
+    _epsilon=1e-6;
+    _indexFacePeriodicSet=false;
+       _nxyz.resize(_spaceDim);
+       _nxyz[0]=nx;
+       _nxyz[1]=ny;
 
-void
-Mesh::setBoundaryNodesFromFaces()
-{
-    for (int iface=0;iface<_boundaryFaceIds.size();iface++)
-    {
-        std::vector< int > nodesID= _faces[_boundaryFaceIds[iface]].getNodesId();
-        int nbNodes = _faces[_boundaryFaceIds[iface]].getNumberOfNodes();
-        for(int inode=0 ; inode<nbNodes ; inode++)
-        {
-            std::vector<int>::const_iterator  it = std::find(_boundaryNodeIds.begin(),_boundaryNodeIds.end(),nodesID[inode]);
-            if( it != _boundaryNodeIds.end() )
-                _boundaryNodeIds.push_back(nodesID[inode]);
-        }
-    }
-}
+       double originPtr[_spaceDim];
+       double dxyzPtr[_spaceDim];
+       mcIdType nodeStrctPtr[_spaceDim];
 
-std::map<int,int>
-Mesh::getIndexFacePeriodic( void ) const
-{
-    return _indexFacePeriodicMap;
-}
+       originPtr[0]=xmin;
+       originPtr[1]=ymin;
+       nodeStrctPtr[0]=nx+1;
+       nodeStrctPtr[1]=ny+1;
+       dxyzPtr[0]=dx;
+       dxyzPtr[1]=dy;
 
-void
-Mesh::setPeriodicFaces(bool check_groups, bool use_central_inversion)
-{
-    if(_indexFacePeriodicSet)
-        return;
-        
-    for (int indexFace=0;indexFace<_boundaryFaceIds.size() ; indexFace++)
-    {
-        Face my_face=_faces[_boundaryFaceIds[indexFace]];
-        int iface_perio=-1;
-        if(_meshDim==1)
-        {
-            for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
-                if(iface!=indexFace)
-                {
-                    iface_perio=_boundaryFaceIds[iface];
-                    break;
-                }
-        }
-        else if(_meshDim==2)
+       _mesh=MEDCouplingIMesh::New(meshName,
+                       _spaceDim,
+                       nodeStrctPtr,
+                       nodeStrctPtr+_spaceDim,
+                       originPtr,
+                       originPtr+_spaceDim,
+                       dxyzPtr,
+                       dxyzPtr+_spaceDim);
+    _meshNotDeleted=true;
+    _isStructured = true;
+
+    if(split_to_triangles_policy==0 || split_to_triangles_policy==1)
         {
-            double x=my_face.x();
-            double y=my_face.y();
-            
-            for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
-            {
-                Face face_i=_faces[_boundaryFaceIds[iface]];
-                double xi=face_i.x();
-                double yi=face_i.y();
-                if (   (abs(y-yi)<_epsilon || abs(x-xi)<_epsilon )// Case of a square geometry
-                    && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked
-                    && ( !use_central_inversion || abs(y+yi) + abs(x+xi)<_epsilon ) // Case of a central inversion
-                    && fabs(my_face.getMeasure() - face_i.getMeasure())<_epsilon
-                    && fabs(my_face.getXN()      + face_i.getXN())<_epsilon
-                    && fabs(my_face.getYN()      + face_i.getYN())<_epsilon
-                    && fabs(my_face.getZN()      + face_i.getZN())<_epsilon )
-                {
-                    iface_perio=_boundaryFaceIds[iface];
-                    break;
-                }
-            }
+            _mesh=_mesh->buildUnstructured();//simplexize is not available for structured meshes
+            _mesh->simplexize(split_to_triangles_policy);
+                       _isStructured = false;
         }
-        else if(_meshDim==3)
+    else if (split_to_triangles_policy != -1)
         {
-            double x=my_face.x();
-            double y=my_face.y();
-            double z=my_face.z();
-        
-            for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
-            {
-                Face face_i=_faces[_boundaryFaceIds[iface]];
-                double xi=face_i.x();
-                double yi=face_i.y();
-                double zi=face_i.z();
-                if ( ((abs(y-yi)<_epsilon && abs(x-xi)<_epsilon) || (abs(x-xi)<_epsilon && abs(z-zi)<_epsilon) || (abs(y-yi)<_epsilon && abs(z-zi)<_epsilon))// Case of a cube geometry
-                    && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked
-                    && ( !use_central_inversion || abs(y+yi) + abs(x+xi) + abs(z+zi)<_epsilon )// Case of a central inversion
-                    && fabs(my_face.getMeasure() - face_i.getMeasure())<_epsilon
-                    && fabs(my_face.getXN()      + face_i.getXN())<_epsilon
-                    && fabs(my_face.getYN()      + face_i.getYN())<_epsilon
-                    && fabs(my_face.getZN()      + face_i.getZN())<_epsilon )
-                {
-                    iface_perio=_boundaryFaceIds[iface];
-                    break;
-                }
-            }  
+            cout<< "split_to_triangles_policy = "<< split_to_triangles_policy << endl;
+            throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : Unknown splitting policy");
         }
-        else
-            throw CdmathException("Mesh::setPeriodicFaces: Mesh dimension should be 1, 2 or 3");
-        
-        if (iface_perio==-1)
-            throw CdmathException("Mesh::setPeriodicFaces: periodic face not found, iface_perio==-1 " );
-        else
-            _indexFacePeriodicMap[_boundaryFaceIds[indexFace]]=iface_perio;
-    }
-    _indexFacePeriodicSet=true;    
+    
+       setMesh();
 }
 
-int
-Mesh::getIndexFacePeriodic(int indexFace, bool check_groups, bool use_central_inversion)
+//----------------------------------------------------------------------
+Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz, int split_to_tetrahedra_policy, std::string meshName)
+//----------------------------------------------------------------------
 {
-       if (!_faces[indexFace].isBorder())
-        {
-            cout<<"Pb with indexFace= "<<indexFace<<endl;
-            throw CdmathException("Mesh::getIndexFacePeriodic: not a border face" );
-        }
-        
-    if(!_indexFacePeriodicSet)
-        setPeriodicFaces(check_groups, use_central_inversion);
-
-    std::map<int,int>::const_iterator  it = _indexFacePeriodicMap.find(indexFace);
-    if( it != _indexFacePeriodicMap.end() )
-        return it->second;
-    else
-    {
-        cout<<"Pb with indexFace= "<<indexFace<<endl;
-        throw CdmathException("Mesh::getIndexFacePeriodic: not a periodic face" );
-    }
-}
+       if(nx<=0 || ny<=0 || nz<=0)
+               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : nx <= 0 or ny <= 0 or nz <= 0");
+       if(xmin>=xmax)
+               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : xmin >= xmax");
+       if(ymin>=ymax)
+               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : ymin >= ymax");
+       if(zmin>=zmax)
+               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : zmin >= zmax");
 
-bool
-Mesh::isBorderNode(int nodeid) const
-{
-       return _nodes[nodeid].isBorder();
-}
+       _spaceDim = 3;
+       _meshDim  = 3;
+    _name=meshName;
+    _epsilon=1e-6;
 
-bool
-Mesh::isBorderFace(int faceid) const
-{
-       return _faces[faceid].isBorder();
-}
+       double dx = (xmax - xmin)/nx ;
+       double dy = (ymax - ymin)/ny ;
+       double dz = (zmax - zmin)/nz ;
 
-std::vector< int > 
-Mesh::getBoundaryFaceIds() const
-{
-    return _boundaryFaceIds;
-}
+       _nxyz.resize(_spaceDim);
+       _nxyz[0]=nx;
+       _nxyz[1]=ny;
+       _nxyz[2]=nz;
 
-std::vector< int > 
-Mesh::getBoundaryNodeIds() const
-{
-    return _boundaryNodeIds;
-}
+       double originPtr[_spaceDim];
+       double dxyzPtr[_spaceDim];
+       mcIdType nodeStrctPtr[_spaceDim];
 
-bool
-Mesh::isTriangular() const
-{
-       return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TRI3;
-}
-bool
-Mesh::isTetrahedral() const
-{
-       return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TETRA4;
-}
-bool
-Mesh::isQuadrangular() const
-{
-       return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_QUAD4;
-}
-bool
-Mesh::isHexahedral() const
-{
-       return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_HEXA8;
-}
-bool
-Mesh::isStructured() const
-{
-       return _isStructured;
-}
+       originPtr[0]=xmin;
+       originPtr[1]=ymin;
+       originPtr[2]=zmin;
+       nodeStrctPtr[0]=nx+1;
+       nodeStrctPtr[1]=ny+1;
+       nodeStrctPtr[2]=nz+1;
+       dxyzPtr[0]=dx;
+       dxyzPtr[1]=dy;
+       dxyzPtr[2]=dz;
 
-std::vector< INTERP_KERNEL::NormalizedCellType > 
-Mesh::getElementTypes() const
-{
-  return _eltsTypes;  
-}
+       _mesh=MEDCouplingIMesh::New(meshName,
+                       _spaceDim,
+                       nodeStrctPtr,
+                       nodeStrctPtr+_spaceDim,
+                       originPtr,
+                       originPtr+_spaceDim,
+                       dxyzPtr,
+                       dxyzPtr+_spaceDim);
+    _meshNotDeleted=true;
+    _isStructured = true;
 
-std::vector< string >
-Mesh::getElementTypesNames() const
-{
-    std::vector< string > result(0);    
-    for(int i=0; i< _eltsTypes.size(); i++)
-    {
-        if( _eltsTypes[i]==INTERP_KERNEL::NORM_POINT1)
-            result.push_back("Points");
-        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_SEG2)
-            result.push_back("Segments");
-        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TRI3)
-            result.push_back("Triangles");
-        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_QUAD4)
-            result.push_back("Quadrangles");
-        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYGON)
-            result.push_back("Polygons");
-        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TETRA4)
-            result.push_back("Tetrahedra");
-        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PYRA5)
-            result.push_back("Pyramids");
-        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PENTA6)
-            result.push_back("Pentahedra");
-        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_HEXA8)
-            result.push_back("Hexahedra");
-        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYHED)
-            result.push_back("Polyhedrons");
-        else
-               {
-                       cout<< "Mesh " + _name + " contains an element of type " <<endl;
-                       cout<< _eltsTypes[i]<<endl;
-                       throw CdmathException("Mesh::getElementTypes : recognised cell med types are NORM_POINT1, NORM_SEG2, NORM_TRI3, NORM_QUAD4, NORM_TETRA4, NORM_PYRA5, NORM_PENTA6, NORM_HEXA8, NORM_POLYGON, NORM_POLYHED");
+    if( split_to_tetrahedra_policy == 0 )
+        {
+            _mesh=_mesh->buildUnstructured();//simplexize is not available for structured meshes
+            _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_5);
+                       _isStructured = false;
         }
-    }
-    return result;
-}
-
-void
-Mesh::setGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh*  mu)
-{
-       int nbCellsSubMesh, nbNodesSubMesh;
-       bool foundFace, foundNode;
-       
-       /* Searching for face groups */
-       vector<string> faceGroups=medmesh->getGroupsNames() ;
-
-       for (unsigned int i=0;i<faceGroups.size();i++ )
-       {
-               string groupName=faceGroups[i];
-               vector<mcIdType> nonEmptyGrp(medmesh->getGrpNonEmptyLevels(groupName));
-               //We check if the group has a relative dimension equal to -1 
-               //before call to the function getGroup(-1,groupName.c_str())
-               vector<mcIdType>::iterator it = find(nonEmptyGrp.begin(), nonEmptyGrp.end(), -1);
-               if (it != nonEmptyGrp.end())
-               {
-                       MEDCouplingUMesh *m=medmesh->getGroup(-1,groupName.c_str());
-            m->unPolyze();
-                       m->sortCellsInMEDFileFrmt( );
-                       nbCellsSubMesh=m->getNumberOfCells();
-                       
-                       _faceGroups.insert(_faceGroups.end(),m);//Vector of group meshes
-                       _faceGroupNames.insert(_faceGroupNames.end(),groupName);//Vector of group names
-                       _faceGroupsIds.insert(_faceGroupsIds.end(),std::vector<int>(nbCellsSubMesh));//Vector of group face Ids. The filling of the face ids takes place below.
-                       
-                       DataArrayDouble *baryCell = m->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;
-                       const double *coorBary=baryCell->getConstPointer();
-                       /* Face identification */
-                       for (int ic(0), k(0); ic<nbCellsSubMesh; ic++, k+=_spaceDim)
-                       {
-                               vector<double> coorBaryXyz(3,0);
-                               for (int d=0; d<_spaceDim; d++)
-                                       coorBaryXyz[d] = coorBary[k+d];
-                               Point p1(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
-
-                               foundFace=false;
-                               for (int iface=0;iface<_numberOfFaces;iface++ )
-                               {
-                                       Point p2=_faces[iface].getBarryCenter();
-                                       if(p1.distance(p2)<_epsilon)
-                                       {
-                                               _faces[iface].setGroupName(groupName);
-                                               _faceGroupsIds[_faceGroupsIds.size()-1][ic]=iface;
-                                               foundFace=true;
-                                               break;
-                                       }
-                               }
-                               if (not foundFace)
-                                       throw CdmathException("No face found for group " + groupName );
-                       }
-                       baryCell->decrRef();
-                       //m->decrRef();
-               }
-       }
-
-       /* Searching for node groups */
-       vector<string> nodeGroups=medmesh->getGroupsOnSpecifiedLev(1) ;
-
-       for (unsigned int i=0;i<nodeGroups.size();i++ )
-       {
-               string groupName=nodeGroups[i];
-               DataArrayIdType * nodeGroup=medmesh->getNodeGroupArr( groupName );
-               const mcIdType *nodeids=nodeGroup->getConstPointer();
-
-               if(nodeids!=NULL)
-               {
-                       _nodeGroups.insert(_nodeGroups.end(),nodeGroup);
-                       _nodeGroupNames.insert(_nodeGroupNames.end(),groupName);
-
-                       nbNodesSubMesh=nodeGroup->getNumberOfTuples();//nodeGroup->getNbOfElems();
-
-                       DataArrayDouble *coo = mu->getCoords() ;
-                       const double *cood=coo->getConstPointer();
-
-                       _nodeGroupsIds.insert(_nodeGroupsIds.end(),std::vector<int>(nbNodesSubMesh));//Vector of boundary faces
-                       /* Node identification */
-                       for (int ic(0); ic<nbNodesSubMesh; ic++)
-                       {
-                               vector<double> coorP(3,0);
-                               for (int d=0; d<_spaceDim; d++)
-                                       coorP[d] = cood[nodeids[ic]*_spaceDim+d];
-                               Point p1(coorP[0],coorP[1],coorP[2]) ;
-
-                               foundNode=false;
-                               for (int inode=0;inode<_numberOfNodes;inode++ )
-                               {
-                                       Point p2=_nodes[inode].getPoint();
-                                       if(p1.distance(p2)<_epsilon)
-                                       {
-                                               _nodes[inode].setGroupName(groupName);
-                                               _nodeGroupsIds[_nodeGroupsIds.size()-1][ic]=inode;
-                                               foundNode=true;
-                                               break;
-                                       }
-                               }
-                               if (not foundNode)
-                                       throw CdmathException("No node found for group " + groupName );
-                       }
-               }
-       }
+    else if( split_to_tetrahedra_policy == 1 )
+        {
+            _mesh=_mesh->buildUnstructured();//simplexize is not available for structured meshes
+            _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_6);
+                       _isStructured = false;
+        }
+    else if ( split_to_tetrahedra_policy != -1 )
+        {
+            cout<< "split_to_tetrahedra_policy = "<< split_to_tetrahedra_policy << endl;
+            throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : splitting policy value should be 0 or 1");
+        }
+    
+       setMesh();
 }
 
 //----------------------------------------------------------------------
@@ -888,13 +496,13 @@ Mesh::setMesh( void )
        const mcIdType *tmpNEI=nodalI->getConstPointer();
 
        _numberOfCells = mu->getNumberOfCells() ;
-       _cells      = new Cell[_numberOfCells] ;
+       _cells      = std::shared_ptr<Cell>(new Cell[_numberOfCells], std::default_delete<Cell[]>()) ;
 
        _numberOfNodes = mu->getNumberOfNodes() ;//This number may include isolated nodes that will not be loaded. The number will be updated during nodes constructions
-       _nodes      = new Node[_numberOfNodes] ;//This array may be resized if isolated nodes are found
+       _nodes      = std::shared_ptr<Node>(new Node[_numberOfNodes], std::default_delete<Node[]>())  ;//This array may be resized if isolated nodes are found
 
        _numberOfFaces = mu2->getNumberOfCells();
-       _faces       = new Face[_numberOfFaces] ;
+       _faces       = std::shared_ptr<Face>(new Face[_numberOfFaces], std::default_delete<Face[]>()) ;
 
     _indexFacePeriodicSet=false;
 
@@ -966,7 +574,7 @@ Mesh::setMesh( void )
                                ci.addFaceId(el,work[el]) ;
                                xn = - xn; yn=-yn; zn=-zn;
                        }
-                       _cells[id] = ci ;
+                       _cells.get()[id] = ci ;
                }
        
                for( int id(0); id<_numberOfFaces; id++ )
@@ -986,7 +594,7 @@ Mesh::setMesh( void )
                        assert( nbCells>0);//To make sure our face is not located on an isolated node
                    
                        Face fi( nbNodes, nbCells, 1.0, p, 1., 0., 0. ) ;
-                       for(int node_id=0; node_id<nbNodes;node_id++)//This loop could b deleted since nbNodes=1. Trying to merge with setMesh
+                       for(int node_id=0; node_id<nbNodes;node_id++)//This loop could be deleted since nbNodes=1. Trying to merge with setMesh
                                fi.addNodeId(node_id,workv[node_id]) ;//global node number
        
                        fi.addCellId(0,workc[0]) ;
@@ -1001,7 +609,7 @@ Mesh::setMesh( void )
                        }
                        if(nbCells==1)
                                _boundaryFaceIds.push_back(id);
-                       _faces[id] = fi ;
+                       _faces.get()[id] = fi ;
                }
        
                int correctNbNodes=0;
@@ -1032,14 +640,18 @@ Mesh::setMesh( void )
                        for( int el=0;el<nbNeighbourNodes;el++ )
                                        vi.addNeighbourNodeId(el,workn[el]) ;//global node number
                                for( int el=0;el<nbFaces;el++ )
-                                       vi.addFaceId(el,workf[el],_faces[workf[el]].isBorder()) ;
+                                       vi.addFaceId(el,workf[el],_faces.get()[workf[el]].isBorder()) ;
                                if(vi.isBorder())
                                        _boundaryNodeIds.push_back(id);
-                               _nodes[id] = vi ;
+                               _nodes.get()[id] = vi ;
                        }
                }
-               if( _numberOfNodes!=correctNbNodes)
+               if( _numberOfNodes!=correctNbNodes)//To do : reduce the size of pointer _nodes
+               {
                        cout<<"Found isolated nodes : correctNbNodes= "<<correctNbNodes<<", _numberOfNodes= "<<_numberOfNodes<<endl;
+                       _numberOfNodes = correctNbNodes;
+                       //memcpy(_nodes,mesh.getNodes(),correctNbNodes*sizeof(*mesh.getNodes())) ;
+               }
        }
        else if(_meshDim==2  || _meshDim==3)
        {
@@ -1160,7 +772,7 @@ Mesh::setMesh( void )
                                        }
                                }
                        }
-                       _cells[id] = ci ;
+                       _cells.get()[id] = ci ;
                }
 
                if(_spaceDim!=_meshDim)
@@ -1219,7 +831,7 @@ Mesh::setMesh( void )
                     }                
             }
             
-                       _faces[id] = fi ;
+                       _faces.get()[id] = fi ;
                }
 
                /*Building mesh nodes, should be done after building mesh faces in order to detect boundary nodes*/
@@ -1264,15 +876,18 @@ Mesh::setMesh( void )
                                        vi.addNeighbourNodeId(el,workn[el]) ;
                                //Detection of border nodes    
                                for( int el=0;el<nbFaces;el++ )
-                                       vi.addFaceId(el,workf[el],_faces[workf[el]].isBorder()) ;
+                                       vi.addFaceId(el,workf[el],_faces.get()[workf[el]].isBorder()) ;
                                if(vi.isBorder())
                                        _boundaryNodeIds.push_back(id);
-                               _nodes[id] = vi ;
+                               _nodes.get()[id] = vi ;
                        }
                }
-               if( _numberOfNodes!=correctNbNodes)
+               if( _numberOfNodes!=correctNbNodes)//To do : reduce the size of pointer _nodes
+               {
                        cout<<"Found isolated nodes : correctNbNodes= "<<correctNbNodes<<", _numberOfNodes= "<<_numberOfNodes<<endl;
-
+                       _numberOfNodes = correctNbNodes;
+               }
+               
                if(_spaceDim==_meshDim)
                        fieldn->decrRef();
                fieldl->decrRef();
@@ -1281,345 +896,568 @@ Mesh::setMesh( void )
        else
                throw CdmathException("Mesh::setMesh space dimension should be 1, 2 or 3");
 
-    //definition of the bounding box for unstructured meshes
-    if(!_isStructured)//Structured case is treated in function readMeshMed
+    //Set boundary groups
+    _faceGroupNames.push_back("Boundary");
+    _nodeGroupNames.push_back("Boundary");
+    _faceGroupsIds.push_back(_boundaryFaceIds);
+    _nodeGroupsIds.push_back(_boundaryNodeIds);
+    if( _meshDim>1 )
+    {    //Set face boundary group
+               _boundaryMesh = mu->computeSkin();
+               _faceGroups.push_back(_boundaryMesh);
+       }
+       else
+               _faceGroups.push_back(NULL);
+    _nodeGroups.push_back(NULL);
+
+    desc->decrRef();
+       descI->decrRef();
+       revDesc->decrRef();
+       revDescI->decrRef();
+       mu2->decrRef();
+       baryCell->decrRef();
+       fields->decrRef();
+       revNode->decrRef();
+       revNodeI->decrRef();
+       revCell->decrRef();
+       revCellI->decrRef();
+
+    if (_meshDim == 3)
+    {
+        revNode2->decrRef();
+        revNodeI2->decrRef();
+        desc2->decrRef();
+        descI2->decrRef();
+        revDesc2->decrRef();
+        revDescI2->decrRef();
+        mu3->decrRef();
+    }
+       
+    return mu;
+}
+
+void
+Mesh::setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup)
+{
+       std::vector< int > faceIds(0);
+       double FX, FY, FZ;
+       Face Fi;
+       
+       /* Construction of the face group */
+       if(isBoundaryGroup)        
+        for(int i=0; i<_boundaryFaceIds.size(); i++)
+        {
+                       Fi=_faces.get()[_boundaryFaceIds[i]];
+                       FX=Fi.x();
+                       FY=Fi.y();
+                       FZ=Fi.z();
+                       if (abs(FX-x)<eps && abs(FY-y)<eps && abs(FZ-z)<eps)
+                       {
+                               faceIds.insert(faceIds.end(),_boundaryFaceIds[i]);
+                               _faces.get()[_boundaryFaceIds[i]].setGroupName(groupName);
+                       }
+        }
+       else
+               for (int iface=0;iface<_numberOfFaces;iface++)
+               {
+                       Fi=_faces.get()[iface];
+                       FX=Fi.x();
+                       FY=Fi.y();
+                       FZ=Fi.z();
+                       if (abs(FX-x)<eps && abs(FY-y)<eps && abs(FZ-z)<eps)
+                       {
+                               faceIds.insert(faceIds.end(),iface);
+                               _faces.get()[iface].setGroupName(groupName);
+                       }
+               }
+
+       if (faceIds.size()>0)
+    {
+               std::vector< std::string >::iterator it = std::find(_faceGroupNames.begin(), _faceGroupNames.end(), groupName);
+               if(it == _faceGroupNames.end())//No group named groupName
+               {
+                       _faceGroupNames.insert(_faceGroupNames.end(),groupName);
+                       _faceGroupsIds.insert(  _faceGroupsIds.end(),faceIds);
+                       _faceGroups.insert(    _faceGroups.end(), NULL);//No mesh created. Create one ?
+               }
+               else
+               {
+                       std::vector< int > faceGroupIds = _faceGroupsIds[it-_faceGroupNames.begin()];
+                       faceGroupIds.insert( faceGroupIds.end(), faceIds.begin(), faceIds.end());
+                       /* Detect and erase duplicates face ids */
+                       sort( faceGroupIds.begin(), faceGroupIds.end() );
+                       faceGroupIds.erase( unique( faceGroupIds.begin(), faceGroupIds.end() ), faceGroupIds.end() );
+                       _faceGroupsIds[it-_faceGroupNames.begin()] = faceGroupIds;
+               }
+       }
+}
+
+void
+Mesh::setGroupAtNodeByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup)
+{
+       std::vector< int > nodeIds(0);
+       double NX, NY, NZ;
+       Node Ni;
+       
+       /* Construction of the node group */
+       if(isBoundaryGroup)        
+        for(int i=0; i<_boundaryNodeIds.size(); i++)
+        {
+                       Ni=_nodes.get()[_boundaryNodeIds[i]];
+                       NX=Ni.x();
+                       NY=Ni.y();
+                       NZ=Ni.z();
+                       if (abs(NX-x)<eps && abs(NY-y)<eps && abs(NZ-z)<eps)
+                       {
+                               nodeIds.insert(nodeIds.end(),_boundaryNodeIds[i]);
+                               _nodes.get()[_boundaryNodeIds[i]].setGroupName(groupName);
+                       }
+        }
+       else
+               for (int inode=0;inode<_numberOfNodes;inode++)
+               {
+                       NX=_nodes.get()[inode].x();
+                       NY=_nodes.get()[inode].y();
+                       NZ=_nodes.get()[inode].z();
+                       if (abs(NX-x)<eps && abs(NY-y)<eps && abs(NZ-z)<eps)
+                       {
+                               nodeIds.insert(nodeIds.end(),inode);
+                               _nodes.get()[inode].setGroupName(groupName);
+                       }
+               }
+
+       if (nodeIds.size()>0)
+    {
+               std::vector< std::string >::iterator it = std::find(_nodeGroupNames.begin(), _nodeGroupNames.end(), groupName);
+               if(it == _nodeGroupNames.end())//No group named groupName
+               {
+                       _nodeGroupNames.insert(_nodeGroupNames.end(),groupName);
+                       _nodeGroupsIds.insert(  _nodeGroupsIds.end(),nodeIds);
+                       _nodeGroups.insert(    _nodeGroups.end(), NULL);//No mesh created. Create one ?
+               }
+               else
+               {
+                       std::vector< int > nodeGroupIds = _nodeGroupsIds[it-_nodeGroupNames.begin()];
+                       nodeGroupIds.insert( nodeGroupIds.end(), nodeIds.begin(), nodeIds.end());
+                       /* Detect and erase duplicates node ids */
+                       sort( nodeGroupIds.begin(), nodeGroupIds.end() );
+                       nodeGroupIds.erase( unique( nodeGroupIds.begin(), nodeGroupIds.end() ), nodeGroupIds.end() );
+                       _nodeGroupsIds[it-_nodeGroupNames.begin()] = nodeGroupIds;
+               }
+    }
+}
+
+void
+Mesh::setGroupAtPlan(double value, int direction, double eps, std::string groupName, bool isBoundaryGroup)
+{
+       std::vector< int > faceIds(0), nodeIds(0);
+       double cord;
+       
+       /* Construction of the face group */    
+       if(isBoundaryGroup)        
+        for(int i=0; i<_boundaryFaceIds.size(); i++)
+        {
+                       cord=_faces.get()[_boundaryFaceIds[i]].getBarryCenter()[direction];
+                       if (abs(cord-value)<eps)
+                       {
+                               faceIds.insert(faceIds.end(),_boundaryFaceIds[i]);
+                               _faces.get()[_boundaryFaceIds[i]].setGroupName(groupName);
+                       }
+        }
+       else
+               for (int iface=0;iface<_numberOfFaces;iface++)
+               {
+                       cord=_faces.get()[iface].getBarryCenter()[direction];
+                       if (abs(cord-value)<eps)
+                       {
+                               faceIds.insert(faceIds.end(),iface);
+                               _faces.get()[iface].setGroupName(groupName);
+                       }
+               }
+
+       /* Construction of the node group */
+       if(isBoundaryGroup)        
+        for(int i=0; i<_boundaryNodeIds.size(); i++)
+        {
+                       cord=_nodes.get()[_boundaryNodeIds[i]].getPoint()[direction];
+                       if (abs(cord-value)<eps)
+                       {
+                               nodeIds.insert(nodeIds.end(),_boundaryNodeIds[i]);
+                               _nodes.get()[_boundaryNodeIds[i]].setGroupName(groupName);
+                       }
+        }
+       else
+               for (int inode=0;inode<_numberOfNodes;inode++)
+               {
+                       cord=_nodes.get()[inode].getPoint()[direction];
+                       if (abs(cord-value)<eps)
+                       {
+                               nodeIds.insert(nodeIds.end(),inode);
+                               _nodes.get()[inode].setGroupName(groupName);
+                       }
+               }
+
+       if (faceIds.size()>0)
+    {
+               std::vector< std::string >::iterator it = std::find(_faceGroupNames.begin(), _faceGroupNames.end(), groupName);
+               if(it == _faceGroupNames.end())//No group named groupName
+               {
+                       _faceGroupNames.insert(_faceGroupNames.end(),groupName);
+                       _faceGroupsIds.insert(  _faceGroupsIds.end(),faceIds);
+                       _faceGroups.insert(    _faceGroups.end(), NULL);//No mesh created. Create one ?
+               }
+               else
+               {cout<<"_faceGroupNames.size()="<<_faceGroupNames.size()<<", _faceGroupsIds.size()="<<_faceGroupsIds.size()<<endl;
+                       std::vector< int > faceGroupIds = _faceGroupsIds[it-_faceGroupNames.begin()];
+                       faceGroupIds.insert( faceGroupIds.end(), faceIds.begin(), faceIds.end());
+                       /* Detect and erase duplicates face ids */
+                       sort( faceGroupIds.begin(), faceGroupIds.end() );
+                       faceGroupIds.erase( unique( faceGroupIds.begin(), faceGroupIds.end() ), faceGroupIds.end() );
+                       _faceGroupsIds[it-_faceGroupNames.begin()] = faceGroupIds;
+               }
+       }
+       if (nodeIds.size()>0)
+    {
+               std::vector< std::string >::iterator it = std::find(_nodeGroupNames.begin(), _nodeGroupNames.end(), groupName);
+               if(it == _nodeGroupNames.end())//No group named groupName
+               {
+                       _nodeGroupNames.insert(_nodeGroupNames.end(),groupName);
+                       _nodeGroupsIds.insert(  _nodeGroupsIds.end(),nodeIds);
+                       _nodeGroups.insert(    _nodeGroups.end(), NULL);//No mesh created. Create one ?
+               }
+               else
+               {
+                       std::vector< int > nodeGroupIds = _nodeGroupsIds[it-_nodeGroupNames.begin()];
+                       nodeGroupIds.insert( nodeGroupIds.end(), nodeIds.begin(), nodeIds.end());
+                       /* Detect and erase duplicates node ids */
+                       sort( nodeGroupIds.begin(), nodeGroupIds.end() );
+                       nodeGroupIds.erase( unique( nodeGroupIds.begin(), nodeGroupIds.end() ), nodeGroupIds.end() );
+                       _nodeGroupsIds[it-_nodeGroupNames.begin()] = nodeGroupIds;
+               }
+    }
+}
+
+void
+Mesh::setBoundaryNodesFromFaces()
+{
+    for (int iface=0;iface<_boundaryFaceIds.size();iface++)
+    {
+        std::vector< int > nodesID= _faces.get()[_boundaryFaceIds[iface]].getNodesId();
+        int nbNodes = _faces.get()[_boundaryFaceIds[iface]].getNumberOfNodes();
+        for(int inode=0 ; inode<nbNodes ; inode++)
+        {
+            std::vector<int>::const_iterator  it = std::find(_boundaryNodeIds.begin(),_boundaryNodeIds.end(),nodesID[inode]);
+            if( it != _boundaryNodeIds.end() )
+                _boundaryNodeIds.push_back(nodesID[inode]);
+        }
+    }
+}
+
+std::map<int,int>
+Mesh::getIndexFacePeriodic( void ) const
+{
+    return _indexFacePeriodicMap;
+}
+
+void
+Mesh::setPeriodicFaces(bool check_groups, bool use_central_inversion)
+{
+    if(_indexFacePeriodicSet)
+        return;
+        
+    for (int indexFace=0;indexFace<_boundaryFaceIds.size() ; indexFace++)
     {
-        double Box0[2*_spaceDim];
-        coo->getMinMaxPerComponent(Box0);
-
-        _xMin=Box0[0];
-        _xMax=Box0[1];
-        if (_spaceDim>=2)
+        Face my_face=_faces.get()[_boundaryFaceIds[indexFace]];
+        int iface_perio=-1;
+        if(_meshDim==1)
         {
-            _yMin=Box0[2];
-            _yMax=Box0[3];
+            for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
+                if(iface!=indexFace)
+                {
+                    iface_perio=_boundaryFaceIds[iface];
+                    break;
+                }
         }
-        if (_spaceDim>=3)
+        else if(_meshDim==2)
         {
-            _zMin=Box0[4];
-            _zMax=Box0[5];
+            double x=my_face.x();
+            double y=my_face.y();
+            
+            for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
+            {
+                Face face_i=_faces.get()[_boundaryFaceIds[iface]];
+                double xi=face_i.x();
+                double yi=face_i.y();
+                if (   (abs(y-yi)<_epsilon || abs(x-xi)<_epsilon )// Case of a square geometry
+                    && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked
+                    && ( !use_central_inversion || abs(y+yi) + abs(x+xi)<_epsilon ) // Case of a central inversion
+                    && fabs(my_face.getMeasure() - face_i.getMeasure())<_epsilon
+                    && fabs(my_face.getXN()      + face_i.getXN())<_epsilon
+                    && fabs(my_face.getYN()      + face_i.getYN())<_epsilon
+                    && fabs(my_face.getZN()      + face_i.getZN())<_epsilon )
+                {
+                    iface_perio=_boundaryFaceIds[iface];
+                    break;
+                }
+            }
         }
+        else if(_meshDim==3)
+        {
+            double x=my_face.x();
+            double y=my_face.y();
+            double z=my_face.z();
+        
+            for (int iface=0;iface<_boundaryFaceIds.size() ; iface++)
+            {
+                Face face_i=_faces.get()[_boundaryFaceIds[iface]];
+                double xi=face_i.x();
+                double yi=face_i.y();
+                double zi=face_i.z();
+                if ( ((abs(y-yi)<_epsilon && abs(x-xi)<_epsilon) || (abs(x-xi)<_epsilon && abs(z-zi)<_epsilon) || (abs(y-yi)<_epsilon && abs(z-zi)<_epsilon))// Case of a cube geometry
+                    && ( !check_groups || my_face.getGroupName()!=face_i.getGroupName()) //In case groups need to be checked
+                    && ( !use_central_inversion || abs(y+yi) + abs(x+xi) + abs(z+zi)<_epsilon )// Case of a central inversion
+                    && fabs(my_face.getMeasure() - face_i.getMeasure())<_epsilon
+                    && fabs(my_face.getXN()      + face_i.getXN())<_epsilon
+                    && fabs(my_face.getYN()      + face_i.getYN())<_epsilon
+                    && fabs(my_face.getZN()      + face_i.getZN())<_epsilon )
+                {
+                    iface_perio=_boundaryFaceIds[iface];
+                    break;
+                }
+            }  
+        }
+        else
+            throw CdmathException("Mesh::setPeriodicFaces: Mesh dimension should be 1, 2 or 3");
+        
+        if (iface_perio==-1)
+            throw CdmathException("Mesh::setPeriodicFaces: periodic face not found, iface_perio==-1 " );
+        else
+            _indexFacePeriodicMap[_boundaryFaceIds[indexFace]]=iface_perio;
     }
-       
-    //Set boundary groups
-    _faceGroupNames.push_back("Boundary");
-    _nodeGroupNames.push_back("Boundary");
-    _faceGroupsIds.push_back(_boundaryFaceIds);
-    _nodeGroupsIds.push_back(_boundaryNodeIds);
-    if( _meshDim>1 )
-    {    //Set face boundary group
-               _boundaryMesh = mu->computeSkin();
-               _faceGroups.push_back(_boundaryMesh);
-       }
-       else
-               _faceGroups.push_back(NULL);
-    _nodeGroups.push_back(NULL);
-
-    desc->decrRef();
-       descI->decrRef();
-       revDesc->decrRef();
-       revDescI->decrRef();
-       mu2->decrRef();
-       baryCell->decrRef();
-       fields->decrRef();
-       revNode->decrRef();
-       revNodeI->decrRef();
-       revCell->decrRef();
-       revCellI->decrRef();
-
-    if (_meshDim == 3)
-    {
-        revNode2->decrRef();
-        revNodeI2->decrRef();
-        desc2->decrRef();
-        descI2->decrRef();
-        revDesc2->decrRef();
-        revDescI2->decrRef();
-        mu3->decrRef();
-    }
-       
-    return mu;
+    _indexFacePeriodicSet=true;    
 }
 
-//----------------------------------------------------------------------
-Mesh::Mesh( std::vector<double> points, std::string meshName )
-//----------------------------------------------------------------------
+int
+Mesh::getIndexFacePeriodic(int indexFace, bool check_groups, bool use_central_inversion)
 {
-    int nx=points.size();
-    
-       if(nx<2)
-               throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : nx < 2, vector should contain at least two values");
-    int i=0;
-    while( i<nx-1 && points[i+1]>points[i] )
-        i++;
-       if( i!=nx-1 )
-    {
-        //cout<< points << endl;
-               throw CdmathException("Mesh::Mesh( vector<double> points, string meshName) : vector values should be sorted");
-    }
-    
-       _spaceDim = 1 ;
-       _meshDim  = 1 ;
-    _name=meshName;
-    _epsilon=1e-6;
-       _xMin=points[0];
-       _xMax=points[nx-1];
-       _yMin=0.;
-       _yMax=0.;
-       _zMin=0.;
-       _zMax=0.;
+       if (!_faces.get()[indexFace].isBorder())
+        {
+            cout<<"Pb with indexFace= "<<indexFace<<endl;
+            throw CdmathException("Mesh::getIndexFacePeriodic: not a border face" );
+        }
+        
+    if(!_indexFacePeriodicSet)
+        setPeriodicFaces(check_groups, use_central_inversion);
 
-    _isStructured = false;
-    _indexFacePeriodicSet=false;
-    
-    MEDCouplingUMesh * mesh1d = MEDCouplingUMesh::New(meshName, 1);
-    mesh1d->allocateCells(nx - 1);
-    double * coords = new double[nx];
-    mcIdType * nodal_con = new mcIdType[2];
-    coords[0]=points[0];
-    for(int i=0; i<nx- 1 ; i++)
+    std::map<int,int>::const_iterator  it = _indexFacePeriodicMap.find(indexFace);
+    if( it != _indexFacePeriodicMap.end() )
+        return it->second;
+    else
     {
-        nodal_con[0]=i;
-        nodal_con[1]=i+1;
-        mesh1d->insertNextCell(INTERP_KERNEL::NORM_SEG2, 2, nodal_con);
-        coords[i+1]=points[i + 1];
+        cout<<"Pb with indexFace= "<<indexFace<<endl;
+        throw CdmathException("Mesh::getIndexFacePeriodic: not a periodic face" );
     }
-    mesh1d->finishInsertingCells();
-
-    DataArrayDouble * coords_arr = DataArrayDouble::New();
-    coords_arr->alloc(nx,1);
-    std::copy(coords,coords+nx,coords_arr->getPointer());
-    mesh1d->setCoords(coords_arr);
-
-    delete [] coords, nodal_con;
-    coords_arr->decrRef();
-
-    _mesh=mesh1d->buildUnstructured();//To enable writeMED. Because we declared the mesh as unstructured, we decide to build the unstructured data (not mandatory)
-    _meshNotDeleted=true;
-
-       setMesh();
 }
 
-//----------------------------------------------------------------------
-Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName )
-//----------------------------------------------------------------------
+bool
+Mesh::isBorderNode(int nodeid) const
 {
-       if(nx<=0)
-               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : nx <= 0");
-       if(xmin>=xmax)
-               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx) : xmin >= xmax");
-
-       double dx = (xmax - xmin)/nx ;
-
-       _spaceDim = 1 ;
-       _meshDim  = 1 ;
-    _name=meshName;
-    _epsilon=1e-6;
-    _isStructured = true;
-    _indexFacePeriodicSet=false;
-
-       _xMin=xmin;
-       _xMax=xmax;
-       _yMin=0.;
-       _yMax=0.;
-       _zMin=0.;
-       _zMax=0.;
-
-       _dxyz.resize(_spaceDim);
-       _dxyz[0]=dx;
-       _nxyz.resize(_spaceDim);
-       _nxyz[0]=nx;
-
-       double *originPtr = new double[_spaceDim];
-       double *dxyzPtr = new double[_spaceDim];
-       mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
-
-       originPtr[0]=xmin;
-       nodeStrctPtr[0]=nx+1;
-       dxyzPtr[0]=dx;
-
-       _mesh=MEDCouplingIMesh::New(meshName,
-                       _spaceDim,
-                       nodeStrctPtr,
-                       nodeStrctPtr+_spaceDim,
-                       originPtr,
-                       originPtr+_spaceDim,
-                       dxyzPtr,
-                       dxyzPtr+_spaceDim);
-    _meshNotDeleted=true;//Because the mesh is structured cartesian : no data in memory. No nodes and cell coordinates stored
-
-       delete [] originPtr;
-       delete [] dxyzPtr;
-       delete [] nodeStrctPtr;
-
-       setMesh();
+       return _nodes.get()[nodeid].isBorder();
 }
 
-//----------------------------------------------------------------------
-Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, int split_to_triangles_policy, std::string meshName)
-//----------------------------------------------------------------------
+bool
+Mesh::isBorderFace(int faceid) const
 {
-       if(nx<=0 || ny<=0)
-               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : nx <= 0 or ny <= 0");
-       if(xmin>=xmax)
-               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : xmin >= xmax");
-       if(ymin>=ymax)
-               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : ymin >= ymax");
-
-       _xMin=xmin;
-       _xMax=xmax;
-       _yMin=ymin;
-       _yMax=ymax;
-       _zMin=0.;
-       _zMax=0.;
-
-
-       double dx = (xmax - xmin)/nx ;
-       double dy = (ymax - ymin)/ny ;
-
-       _spaceDim = 2 ;
-       _meshDim  = 2 ;
-    _name=meshName;
-    _epsilon=1e-6;
-    _indexFacePeriodicSet=false;
-       _nxyz.resize(_spaceDim);
-       _nxyz[0]=nx;
-       _nxyz[1]=ny;
+       return _faces.get()[faceid].isBorder();
+}
 
-       _dxyz.resize(_spaceDim);
-       _dxyz[0]=dx;
-       _dxyz[1]=dy;
+std::vector< int > 
+Mesh::getBoundaryFaceIds() const
+{
+    return _boundaryFaceIds;
+}
 
-       double *originPtr = new double[_spaceDim];
-       double *dxyzPtr   = new double[_spaceDim];
-       mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
+std::vector< int > 
+Mesh::getBoundaryNodeIds() const
+{
+    return _boundaryNodeIds;
+}
 
-       originPtr[0]=xmin;
-       originPtr[1]=ymin;
-       nodeStrctPtr[0]=nx+1;
-       nodeStrctPtr[1]=ny+1;
-       dxyzPtr[0]=dx;
-       dxyzPtr[1]=dy;
+bool
+Mesh::isTriangular() const
+{
+       return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TRI3;
+}
+bool
+Mesh::isTetrahedral() const
+{
+       return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_TETRA4;
+}
+bool
+Mesh::isQuadrangular() const
+{
+       return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_QUAD4;
+}
+bool
+Mesh::isHexahedral() const
+{
+       return _eltsTypes.size()==1 && _eltsTypes[0]==INTERP_KERNEL::NORM_HEXA8;
+}
+bool
+Mesh::isStructured() const
+{
+       return _isStructured;
+}
 
-       _mesh=MEDCouplingIMesh::New(meshName,
-                       _spaceDim,
-                       nodeStrctPtr,
-                       nodeStrctPtr+_spaceDim,
-                       originPtr,
-                       originPtr+_spaceDim,
-                       dxyzPtr,
-                       dxyzPtr+_spaceDim);
-    _meshNotDeleted=true;//Because the mesh is structured cartesian : no data in memory. No nodes and cell coordinates stored
-    _isStructured = true;
+std::vector< INTERP_KERNEL::NormalizedCellType > 
+Mesh::getElementTypes() const
+{
+  return _eltsTypes;  
+}
 
-    if(split_to_triangles_policy==0 || split_to_triangles_policy==1)
-        {
-            _mesh=_mesh->buildUnstructured();
-            _mesh->simplexize(split_to_triangles_policy);
-            _meshNotDeleted=true;//Now the mesh is unstructured and stored with nodes and cell coordinates
-                       _isStructured = false;
-        }
-    else if (split_to_triangles_policy != -1)
-        {
-            cout<< "split_to_triangles_policy = "<< split_to_triangles_policy << endl;
-            throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny) : Unknown splitting policy");
+std::vector< string >
+Mesh::getElementTypesNames() const
+{
+    std::vector< string > result(0);    
+    for(int i=0; i< _eltsTypes.size(); i++)
+    {
+        if( _eltsTypes[i]==INTERP_KERNEL::NORM_POINT1)
+            result.push_back("Points");
+        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_SEG2)
+            result.push_back("Segments");
+        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TRI3)
+            result.push_back("Triangles");
+        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_QUAD4)
+            result.push_back("Quadrangles");
+        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYGON)
+            result.push_back("Polygons");
+        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_TETRA4)
+            result.push_back("Tetrahedra");
+        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PYRA5)
+            result.push_back("Pyramids");
+        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_PENTA6)
+            result.push_back("Pentahedra");
+        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_HEXA8)
+            result.push_back("Hexahedra");
+        else if( _eltsTypes[i]==INTERP_KERNEL::NORM_POLYHED)
+            result.push_back("Polyhedrons");
+        else
+               {
+                       cout<< "Mesh " + _name + " contains an element of type " <<endl;
+                       cout<< _eltsTypes[i]<<endl;
+                       throw CdmathException("Mesh::getElementTypes : recognised cell med types are NORM_POINT1, NORM_SEG2, NORM_TRI3, NORM_QUAD4, NORM_TETRA4, NORM_PYRA5, NORM_PENTA6, NORM_HEXA8, NORM_POLYGON, NORM_POLYHED");
         }
-        
-       delete [] originPtr;
-       delete [] dxyzPtr;
-       delete [] nodeStrctPtr;
-    
-       setMesh();
+    }
+    return result;
 }
 
-//----------------------------------------------------------------------
-Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz, int split_to_tetrahedra_policy, std::string meshName)
-//----------------------------------------------------------------------
+void
+Mesh::setFaceGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh*  mu)
 {
-       if(nx<=0 || ny<=0 || nz<=0)
-               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : nx <= 0 or ny <= 0 or nz <= 0");
-       if(xmin>=xmax)
-               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : xmin >= xmax");
-       if(ymin>=ymax)
-               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : ymin >= ymax");
-       if(zmin>=zmax)
-               throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : zmin >= zmax");
-
-       _spaceDim = 3;
-       _meshDim  = 3;
-    _name=meshName;
-    _epsilon=1e-6;
-       _xMin=xmin;
-       _xMax=xmax;
-       _yMin=ymin;
-       _yMax=ymax;
-       _zMin=zmin;
-       _zMax=zmax;
+       int nbCellsSubMesh;
+       bool foundFace;
+       
+       /* Searching for face groups */
+       vector<string> faceGroups=medmesh->getGroupsNames() ;
 
-       double dx = (xmax - xmin)/nx ;
-       double dy = (ymax - ymin)/ny ;
-       double dz = (zmax - zmin)/nz ;
+       for (unsigned int i=0;i<faceGroups.size();i++ )
+       {
+               string groupName=faceGroups[i];
+               vector<mcIdType> nonEmptyGrp(medmesh->getGrpNonEmptyLevels(groupName));
+               //We check if the group has a relative dimension equal to -1 
+               //before call to the function getGroup(-1,groupName.c_str())
+               vector<mcIdType>::iterator it = find(nonEmptyGrp.begin(), nonEmptyGrp.end(), -1);
+               if (it != nonEmptyGrp.end())
+               {
+                       MEDCouplingUMesh *m=medmesh->getGroup(-1,groupName.c_str());
+            m->unPolyze();
+                       nbCellsSubMesh=m->getNumberOfCells();
+                       
+                       _faceGroups.insert(_faceGroups.end(),m);//Vector of group meshes
+                       _faceGroupNames.insert(_faceGroupNames.end(),groupName);//Vector of group names
+                       _faceGroupsIds.insert(_faceGroupsIds.end(),std::vector<int>(nbCellsSubMesh));//Vector of group face Ids. The filling of the face ids takes place below.
+                       
+                       DataArrayDouble *baryCell = m->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;
+                       const double *coorBary=baryCell->getConstPointer();
+                       // Face identification 
+                       for (int ic(0), k(0); ic<nbCellsSubMesh; ic++, k+=_spaceDim)
+                       {
+                               vector<double> coorBaryXyz(3,0);
+                               for (int d=0; d<_spaceDim; d++)
+                                       coorBaryXyz[d] = coorBary[k+d];
+                               Point p1(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
 
-       _dxyz.resize(_spaceDim);
-       _dxyz[0]=dx;
-       _dxyz[1]=dy;
-       _dxyz[2]=dz;
+                               foundFace=false;
+                               for (int iface=0;iface<_numberOfFaces;iface++ )
+                               {
+                                       Point p2=_faces.get()[iface].getBarryCenter();
+                                       if(p1.distance(p2)<_epsilon)
+                                       {
+                                               _faces.get()[iface].setGroupName(groupName);
+                                               _faceGroupsIds[_faceGroupsIds.size()-1][ic]=iface;
+                                               foundFace=true;
+                                               break;
+                                       }
+                               }
+                               if (not foundFace)
+                                       throw CdmathException("No face found for group " + groupName );
+                       }
+                       baryCell->decrRef();
+                       //m->decrRef();
+               }
+       }
+}
+void
+Mesh::setNodeGroups( const MEDFileMesh* medmesh, MEDCouplingUMesh*  mu)
+{
+       int  nbNodesSubMesh;
+       bool foundNode;
+       
+       /* Searching for node groups */
+       vector<string> nodeGroups=medmesh->getGroupsOnSpecifiedLev(1) ;
 
-       _nxyz.resize(_spaceDim);
-       _nxyz[0]=nx;
-       _nxyz[1]=ny;
-       _nxyz[2]=nz;
+       for (unsigned int i=0;i<nodeGroups.size();i++ )
+       {
+               string groupName=nodeGroups[i];
+               DataArrayIdType * nodeGroup=medmesh->getNodeGroupArr( groupName );
+               const mcIdType *nodeids=nodeGroup->getConstPointer();
 
-       double *originPtr = new double[_spaceDim];
-       double *dxyzPtr = new double[_spaceDim];
-       mcIdType *nodeStrctPtr = new mcIdType[_spaceDim];
+               if(nodeids!=NULL)
+               {
+                       _nodeGroups.insert(_nodeGroups.end(),nodeGroup);
+                       _nodeGroupNames.insert(_nodeGroupNames.end(),groupName);
 
-       originPtr[0]=xmin;
-       originPtr[1]=ymin;
-       originPtr[2]=zmin;
-       nodeStrctPtr[0]=nx+1;
-       nodeStrctPtr[1]=ny+1;
-       nodeStrctPtr[2]=nz+1;
-       dxyzPtr[0]=dx;
-       dxyzPtr[1]=dy;
-       dxyzPtr[2]=dz;
+                       nbNodesSubMesh=nodeGroup->getNumberOfTuples();//nodeGroup->getNbOfElems();
 
-       _mesh=MEDCouplingIMesh::New(meshName,
-                       _spaceDim,
-                       nodeStrctPtr,
-                       nodeStrctPtr+_spaceDim,
-                       originPtr,
-                       originPtr+_spaceDim,
-                       dxyzPtr,
-                       dxyzPtr+_spaceDim);
-    _meshNotDeleted=true;//Because the mesh is structured cartesian : no data in memory. Nno nodes and cell coordinates stored
-    _isStructured = true;
+                       DataArrayDouble *coo = mu->getCoords() ;
+                       const double *cood=coo->getConstPointer();
 
-    if( split_to_tetrahedra_policy == 0 )
-        {
-            _mesh=_mesh->buildUnstructured();
-            _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_5);
-            _meshNotDeleted=true;//Now the mesh is unstructured and stored with nodes and cell coordinates
-                       _isStructured = false;
-        }
-    else if( split_to_tetrahedra_policy == 1 )
-        {
-            _mesh=_mesh->buildUnstructured();
-            _mesh->simplexize(INTERP_KERNEL::PLANAR_FACE_6);
-            _meshNotDeleted=true;//Now the mesh is unstructured and stored with nodes and cell coordinates
-                       _isStructured = false;
-        }
-    else if ( split_to_tetrahedra_policy != -1 )
-        {
-            cout<< "split_to_tetrahedra_policy = "<< split_to_tetrahedra_policy << endl;
-            throw CdmathException("Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, double zmin, double zmax, int nz) : splitting policy value should be 0 or 1");
-        }
+                       _nodeGroupsIds.insert(_nodeGroupsIds.end(),std::vector<int>(nbNodesSubMesh));//Vector of boundary faces
+                       /* Node identification */
+                       for (int ic(0); ic<nbNodesSubMesh; ic++)
+                       {
+                               vector<double> coorP(3,0);
+                               for (int d=0; d<_spaceDim; d++)
+                                       coorP[d] = cood[nodeids[ic]*_spaceDim+d];
+                               Point p1(coorP[0],coorP[1],coorP[2]) ;
 
-       delete [] originPtr;
-       delete [] dxyzPtr;
-       delete [] nodeStrctPtr;
-    
-       setMesh();
+                               foundNode=false;
+                               for (int inode=0;inode<_numberOfNodes;inode++ )
+                               {
+                                       Point p2=_nodes.get()[inode].getPoint();
+                                       if(p1.distance(p2)<_epsilon)
+                                       {
+                                               _nodes.get()[inode].setGroupName(groupName);
+                                               _nodeGroupsIds[_nodeGroupsIds.size()-1][ic]=inode;
+                                               foundNode=true;
+                                               break;
+                                       }
+                               }
+                               if (not foundNode)
+                                       throw CdmathException("No node found for group " + groupName );
+                       }
+               }
+       }
 }
 
 //----------------------------------------------------------------------
@@ -1638,15 +1476,6 @@ Mesh::getMeshDimension( void )  const
        return _meshDim ;
 }
 
-std::vector<double>
-Mesh::getDXYZ() const
-{
-    if(!_isStructured)
-               throw CdmathException("std::vector<double> Mesh::getDXYZ() : dx,dy and dz are defined only for structured meshes !");
-
-       return _dxyz;
-}
-
 std::vector<mcIdType>
 Mesh::getCellGridStructure() const
 {
@@ -1674,8 +1503,8 @@ Mesh::getNy( void )  const
 {
     if(!_isStructured)
                throw CdmathException("int Mesh::getNy( void ) : Ny is defined only for structured meshes !");
-       if(_spaceDim < 2)
-               throw CdmathException("int double& Field::operator[ielem] : Ny is not defined in dimension < 2!");
+       if(_meshDim < 2)
+               throw CdmathException("int Mesh::getNy( void ) : Ny is not defined for mesh dimension < 2!");
 
        return _nxyz[1];
 }
@@ -1687,8 +1516,8 @@ Mesh::getNz( void )  const
 {
     if(!_isStructured)
                throw CdmathException("int Mesh::getNz( void ) : Nz is defined only for structured meshes !");
-       if(_spaceDim < 3)
-               throw CdmathException("int Mesh::getNz( void ) : Nz is not defined in dimension < 3!");
+       if(_meshDim < 3)
+               throw CdmathException("int Mesh::getNz( void ) : Nz is not defined for mesh dimension < 3!");
 
        return _nxyz[2];
 }
@@ -1697,8 +1526,11 @@ Mesh::getNz( void )  const
 double
 Mesh::getXMin( void )  const
 //----------------------------------------------------------------------
-{        
-       return _xMin ;
+{
+       double Box0[2*_meshDim];
+    _mesh->getBoundingBox(Box0);
+
+       return Box0[0] ;
 }
 
 //----------------------------------------------------------------------
@@ -1706,7 +1538,10 @@ double
 Mesh::getXMax( void )  const
 //----------------------------------------------------------------------
 {
-       return _xMax ;
+       double Box0[2*_meshDim];
+    _mesh->getBoundingBox(Box0);
+
+       return Box0[1] ;
 }
 
 //----------------------------------------------------------------------
@@ -1714,7 +1549,13 @@ double
 Mesh::getYMin( void )  const
 //----------------------------------------------------------------------
 {
-       return _yMin ;
+       if(_meshDim<2)
+               throw CdmathException("Mesh::getYMin : dimension should be >=2");
+               
+       double Box0[2*_meshDim];
+    _mesh->getBoundingBox(Box0);
+
+       return Box0[2] ;
 }
 
 //----------------------------------------------------------------------
@@ -1722,7 +1563,13 @@ double
 Mesh::getYMax( void )  const
 //----------------------------------------------------------------------
 {
-       return _yMax ;
+       if(_meshDim<2)
+               throw CdmathException("Mesh::getYMax : dimension should be >=2");
+               
+       double Box0[2*_meshDim];
+    _mesh->getBoundingBox(Box0);
+
+       return Box0[3] ;
 }
 
 //----------------------------------------------------------------------
@@ -1730,7 +1577,13 @@ double
 Mesh::getZMin( void )  const
 //----------------------------------------------------------------------
 {
-       return _zMin ;
+       if(_meshDim<3)
+               throw CdmathException("Mesh::getZMin : dimension should be 3");
+               
+       double Box0[2*_meshDim];
+    _mesh->getBoundingBox(Box0);
+
+       return Box0[4] ;
 }
 
 //----------------------------------------------------------------------
@@ -1738,7 +1591,13 @@ double
 Mesh::getZMax( void )  const
 //----------------------------------------------------------------------
 {
-       return _zMax ;
+       if(_meshDim<3)
+               throw CdmathException("Mesh::getZMax : dimension should be 3");
+
+       double Box0[2*_meshDim];
+    _mesh->getBoundingBox(Box0);
+
+       return Box0[5] ;
 }
 
 //----------------------------------------------------------------------
@@ -1782,7 +1641,7 @@ Mesh::getNumberOfEdges ( void ) const
 }
 
 //----------------------------------------------------------------------
-Face*
+std::shared_ptr<Face>
 Mesh::getFaces ( void )  const
 //----------------------------------------------------------------------
 {
@@ -1790,7 +1649,7 @@ Mesh::getFaces ( void )  const
 }
 
 //----------------------------------------------------------------------
-Cell*
+std::shared_ptr<Cell>
 Mesh::getCells ( void ) const
 //----------------------------------------------------------------------
 {
@@ -1802,7 +1661,7 @@ Cell&
 Mesh::getCell ( int i ) const
 //----------------------------------------------------------------------
 {
-       return _cells[i] ;
+       return _cells.get()[i] ;
 }
 
 //----------------------------------------------------------------------
@@ -1810,7 +1669,7 @@ Face&
 Mesh::getFace ( int i ) const
 //----------------------------------------------------------------------
 {
-       return _faces[i] ;
+       return _faces.get()[i] ;
 }
 
 //----------------------------------------------------------------------
@@ -1818,11 +1677,11 @@ Node&
 Mesh::getNode ( int i ) const
 //----------------------------------------------------------------------
 {
-       return _nodes[i] ;
+       return _nodes.get()[i] ;
 }
 
 //----------------------------------------------------------------------
-Node*
+std::shared_ptr<Node>
 Mesh::getNodes ( void )  const
 //----------------------------------------------------------------------
 {
@@ -1871,49 +1730,16 @@ Mesh::operator= ( const Mesh& mesh )
     _meshNotDeleted = mesh.meshNotDeleted();
     
     if(_isStructured)
-    {
         _nxyz = mesh.getCellGridStructure() ;
-        _dxyz=mesh.getDXYZ();
-        _xMin=mesh.getXMin();
-        _xMax=mesh.getXMax();
-        _yMin=mesh.getYMin();
-        _yMax=mesh.getYMax();
-        _zMin=mesh.getZMin();
-        _zMax=mesh.getZMax();
-    }
+
        _faceGroupNames = mesh.getNameOfFaceGroups() ;
        _faceGroups = mesh.getFaceGroups() ;
        _nodeGroupNames = mesh.getNameOfNodeGroups() ;
        _nodeGroups = mesh.getNodeGroups() ;
 
-       if (_nodes)
-       {
-               delete [] _nodes ;
-               _nodes=NULL;
-       }
-       if (_faces)
-       {
-               delete [] _faces ;
-               _faces=NULL;
-       }
-       if (_cells)
-       {
-               delete [] _cells ;
-               _cells=NULL;
-       }
-
-       _nodes   = new Node[_numberOfNodes] ;
-       _faces   = new Face[_numberOfFaces] ;
-       _cells   = new Cell[_numberOfCells] ;
-
-       for (int i=0;i<_numberOfNodes;i++)
-               _nodes[i]=mesh.getNode(i);
-
-       for (int i=0;i<_numberOfFaces;i++)
-               _faces[i]=mesh.getFace(i);
-
-       for (int i=0;i<_numberOfCells;i++)
-               _cells[i]=mesh.getCell(i);
+       _nodes   = mesh.getNodes() ;
+       _faces   = mesh.getFaces() ;
+       _cells   = mesh.getCells() ;
 
     _indexFacePeriodicSet= mesh.isIndexFacePeriodicSet();
     if(_indexFacePeriodicSet)
@@ -1925,8 +1751,8 @@ Mesh::operator= ( const Mesh& mesh )
     _eltsTypes=mesh.getElementTypes();
     _eltsTypesNames=mesh.getElementTypesNames();
 
-       MCAuto<MEDCouplingMesh> m1=mesh.getMEDCouplingMesh()->deepCopy();
-       _mesh=m1;
+       _mesh=mesh.getMEDCouplingMesh()->clone(false);//Clone because you will need to buildUnstructured. No deep copy : it is assumed node coordinates and cell connectivity will not change
+
        return *this;
 }
 
@@ -2008,11 +1834,11 @@ Mesh::getMaxNbNeighbours(EntityType type) const
                int nbNeib;//local number of neighbours
         for(int i=0; i<_numberOfCells; i++)
         {
-            Cell Ci = _cells[i];
+            Cell Ci = _cells.get()[i];
             //Careful with mesh with junctions : do not just take Ci.getNumberOfFaces()
             nbNeib=0;
             for(int j=0; j<Ci.getNumberOfFaces(); j++)
-                nbNeib+=_faces[Ci.getFacesId()[j]].getNumberOfCells()-1;//Without junction this would be +=1
+                nbNeib+=_faces.get()[Ci.getFacesId()[j]].getNumberOfCells()-1;//Without junction this would be +=1
             
             if(result < nbNeib)
                 result=nbNeib;
@@ -2021,8 +1847,8 @@ Mesh::getMaxNbNeighbours(EntityType type) const
     else if(type==NODES)
        {
         for(int i=0; i<_numberOfNodes; i++)
-            if(result < _nodes[i].getNumberOfEdges())
-                result=_nodes[i].getNumberOfEdges();
+            if(result < _nodes.get()[i].getNumberOfEdges())
+                result=_nodes.get()[i].getNumberOfEdges();
        }
     else
                throw CdmathException("Mesh::getMaxNbNeighbours : entity type is not accepted. Should be CELLS or NODES");
@@ -2034,8 +1860,8 @@ void
 Mesh::writeVTK ( const std::string fileName ) const
 //----------------------------------------------------------------------
 {
-       if( !_isStructured && !_meshNotDeleted )
-               throw CdmathException("Mesh::writeVTK : Cannot save mesh : no MEDCouplingUMesh loaded");
+       if( !_meshNotDeleted )
+               throw CdmathException("Mesh::writeVTK : Cannot save mesh : no MEDCouplingMesh loaded (may be deleted)");
                
        string fname=fileName+".vtu";
        _mesh->writeVTK(fname.c_str()) ;
@@ -2043,15 +1869,23 @@ Mesh::writeVTK ( const std::string fileName ) const
 
 //----------------------------------------------------------------------
 void
-Mesh::writeMED ( const std::string fileName ) const
+Mesh::writeMED ( const std::string fileName, bool fromScratch ) const
 //----------------------------------------------------------------------
 {
        if( !_meshNotDeleted )
-               throw CdmathException("Mesh::writeMED : Cannot save mesh : no MEDCouplingUMesh loaded");
-               
+               throw CdmathException("Mesh::writeMED : Cannot save mesh : no MEDCouplingMesh loaded (may be deleted)");
+
        string fname=fileName+".med";
-       MEDCouplingUMesh* mu=_mesh->buildUnstructured();
-       MEDCoupling::WriteUMesh(fname.c_str(),mu,true);
+       if(_isStructured)//Check if we have a medcouplingimesh that can't be written directly
+       {
+               const MEDCoupling::MEDCouplingIMesh* iMesh = dynamic_cast< const MEDCoupling::MEDCouplingIMesh* > ((const MEDCoupling::MEDCouplingMesh*) _mesh);
+               if(iMesh)//medcouplingimesh : Use convertToCartesian in order to write mesh
+                       MEDCoupling::WriteMesh(fname.c_str(),iMesh->convertToCartesian(), fromScratch);
+               else//medcouplingcmesh : save directly
+                       MEDCoupling::WriteMesh(fname.c_str(),_mesh, fromScratch);
+       }
+       else
+               MEDCoupling::WriteMesh(fname.c_str(),_mesh, fromScratch);
 
        /* Try to save mesh groups */
        //MEDFileUMesh meshMEDFile;
@@ -2110,7 +1944,7 @@ Mesh::setNodeGroupByIds(std::vector< int > nodeIds, std::string groupName)
         getNode(nodeIds[i]).setGroupName(groupName);
 }
 
-void Mesh::deleteMEDCouplingUMesh()
+void Mesh::deleteMEDCouplingMesh()
 { 
        if(_meshNotDeleted) 
        {