]> SALOME platform Git repositories - tools/solverlab.git/commitdiff
Salome HOME
Improved Mesh class, especially when meshDim=1 (graphs)
authormichael <michael@is242589.intra.cea.fr>
Fri, 10 Dec 2021 13:18:22 +0000 (14:18 +0100)
committermichael <michael@is242589.intra.cea.fr>
Fri, 10 Dec 2021 13:18:22 +0000 (14:18 +0100)
CDMATH/mesh/inc/Mesh.hxx
CDMATH/mesh/src/Mesh.cxx

index a7ea9c9471064f8e2f6b608e1e6ac25bd98daa37..99ce471b2ace22790187b078a8591e02e7de982a 100644 (file)
@@ -106,7 +106,7 @@ public: //----------------------------------------------------------------
         * @param nx : Number of cells in x direction
         * @param ny : Number of cells in y direction
         * @param nz : Number of cells in z direction
-     * @param split_to_tetrahedra_policy : each cuboid will be split into 5 tetrahedra if value is INTERP_KERNEL::PLANAR_FACE_5 or 6 tetrahedra if the value is INTERP_KERNEL::PLANAR_FACE_6
+     * @param split_to_tetrahedra_policy : each cuboid will be split into 5 tetrahedra if value is 0 or 6 tetrahedra if the value is 1
         */
        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") ;
 
@@ -256,10 +256,15 @@ public: //----------------------------------------------------------------
        MEDCoupling::MCAuto<MEDCoupling::MEDCouplingMesh> getMEDCouplingMesh ( void )  const ;
 
        /**
-        * \brief computes the skin surrounding the mesh
+        * \brief return the skin surrounding the mesh
         */
        Mesh getBoundaryMesh ( void )  const ;
 
+       /**
+        * \brief return a group surrounding the mesh
+        */
+       Mesh getBoundaryGroupMesh ( std::string groupName, int nth_group_match = 0 )  const ;
+
        /**
         * \brief return the list of face group names
         * return _faceGroupNames
@@ -303,11 +308,11 @@ public: //----------------------------------------------------------------
      /** 
       * @return list of face group Ids
       */
-    std::vector< int > getGroupFaceIds(std::string groupName) const;
+    std::vector< int > getFaceGroupIds(std::string groupName, bool isBoundaryGroup=true) const;
     /**
      * @return list of node group Ids
      * */
-    std::vector< int > getGroupNodeIds(std::string groupName) const;
+    std::vector< int > getNodeGroupIds(std::string groupName, bool isBoundaryGroup=true) const;
  
        /**
         * \brief write mesh in the VTK format
@@ -319,9 +324,11 @@ public: //----------------------------------------------------------------
         */
        void writeMED ( const std::string fileName ) const ;
 
-       void setGroupAtPlan(double value, int direction, double eps, std::string groupName) ;
+       void setGroupAtPlan(double value, int direction, double eps, std::string groupName, bool isBoundaryGroup=true) ;
 
-       void setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName) ;
+       void setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup=true) ;
+
+       void setGroupAtNodeByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup=true) ;
 
        void setFaceGroupByIds(std::vector< int > faceIds, std::string groupName) ;
 
@@ -349,9 +356,9 @@ public: //----------------------------------------------------------------
        double getComparisonEpsilon() const {return _epsilon;};
        void setComparisonEpsilon(double epsilon){ _epsilon=epsilon;};
     // Quick comparison of two meshes to see if they are identical with high probability (three cells are compared)
-    void checkFastEquivalWith( Mesh m) const { return getMEDCouplingMesh()->checkFastEquivalWith(m.getMEDCouplingMesh(),1e-6);};
+    void checkFastEquivalWith( Mesh m) const { return getMEDCouplingMesh()->checkFastEquivalWith(m.getMEDCouplingMesh(),_epsilon);};
     // Deep comparison of two meshes to see if they are identical Except for their names and units
-    bool isEqualWithoutConsideringStr( Mesh m) const { return getMEDCouplingMesh()->isEqualWithoutConsideringStr(m.getMEDCouplingMesh(),1e-6);};
+    bool isEqualWithoutConsideringStr( Mesh m) const { return getMEDCouplingMesh()->isEqualWithoutConsideringStr(m.getMEDCouplingMesh(),_epsilon);};
 
     std::vector< std::string > getElementTypesNames() const ;
        /**
@@ -364,12 +371,29 @@ public: //----------------------------------------------------------------
         */
     int getMaxNbNeighbours(EntityType type) const;
     
+    /** 
+     * \brief Delete the medcoupling mesh to save memory space
+     */
+    void deleteMEDCouplingUMesh();
+    
+    /** 
+     * \brief Returns true iff an unstructured mesh has been loaded
+     */
+     bool meshNotDeleted() const {return _meshNotDeleted;}
+    
 private: //----------------------------------------------------------------
 
        MEDCoupling::MEDCouplingUMesh*  setMesh( void ) ;
+       void setGroups( const MEDCoupling::MEDFileUMesh* medmesh, MEDCoupling::MEDCouplingUMesh*  mu) ;//Read all face and node group
+       void addNewFaceGroup( const MEDCoupling::MEDCouplingUMesh *m);//adds one face group in the vectors _faceGroups, _faceGroupNames and _faceGroupIds
+       
+       /*
+        * The MEDCoupling mesh
+        */
+       MEDCoupling::MCAuto<MEDCoupling::MEDCouplingMesh> _mesh;
 
-       void setGroups( const MEDCoupling::MEDFileUMesh* medmesh, MEDCoupling::MEDCouplingUMesh*  mu) ;
-
+       bool _meshNotDeleted;
+       
     std::string _name;
     
        /**
@@ -456,11 +480,20 @@ private: //----------------------------------------------------------------
         * The list of node groups.
         */
        std::vector<MEDCoupling::DataArrayIdType *> _nodeGroups;
+       
        /*
-        * The mesh MEDCoupling
+        * The list of face id in each face groups.
+        */
+       std::vector< std::vector<int> > _faceGroupsIds;
+       
+       /*
+        * The list of node id in each node groups.
+        */
+       std::vector< std::vector<int> > _nodeGroupsIds;
+       
+       /*
+        * Elements types (SEG2, TRI3, QUAD4, HEXA6 ...)
         */
-       MEDCoupling::MCAuto<MEDCoupling::MEDCouplingMesh> _mesh;
-
        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;    
@@ -476,7 +509,7 @@ private: //----------------------------------------------------------------
     /* List of boundary nodes*/
     std::vector< int > _boundaryNodeIds;
     /* Boundary mesh */
-    const MEDCoupling::MEDCouplingUMesh * _boundaryMesh;
+    MEDCoupling::MEDCouplingUMesh * _boundaryMesh;
     
     double _epsilon;
 };
index 1b1d4fe8722aec4592cc647f3a2bcc8b6f7d4fd8..44ca0405d8b8c3bd2f524ca884452a2915bc3040 100644 (file)
@@ -20,6 +20,7 @@
 #include <iterator>
 #include <algorithm> 
 #include <string.h>
+#include <assert.h>
 
 using namespace MEDCoupling;
 using namespace std;
@@ -29,6 +30,7 @@ Mesh::Mesh( void )
 //----------------------------------------------------------------------
 {
        _mesh=NULL;
+    _meshNotDeleted=false;
        _cells=NULL;
        _nodes=NULL;
        _faces=NULL;
@@ -47,12 +49,18 @@ Mesh::Mesh( void )
        _zMax=0.;
     _nxyz.resize(0);
     _dxyz.resize(0.);
+       _boundaryMesh=NULL;
+    _boundaryFaceIds.resize(0);
+    _boundaryNodeIds.resize(0);
        _faceGroupNames.resize(0);
        _faceGroups.resize(0);
+       _faceGroupsIds.resize(0);
        _nodeGroupNames.resize(0);
        _nodeGroups.resize(0);
+       _nodeGroupsIds.resize(0);
     _indexFacePeriodicSet=false;
     _name="";
+    _epsilon=1e-6;
 }
 
 //----------------------------------------------------------------------
@@ -62,9 +70,12 @@ Mesh::~Mesh( void )
        delete [] _cells;
        delete [] _nodes;
        delete [] _faces;
+       
        //for(int i=0; i< _faceGroups.size(); i++)
        //      _faceGroups[i]->decrRef();
        //      _nodeGroups[i]->decrRef();
+       if( _meshNotDeleted)
+               (_mesh.retn())->decrRef();
 }
 
 std::string 
@@ -77,12 +88,12 @@ Mesh::Mesh( const MEDCoupling::MEDCouplingIMesh* mesh )
 {
        _spaceDim=mesh->getSpaceDimension();
        _meshDim=mesh->getMeshDimension();
-    _isStructured=true;
        _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];
@@ -116,10 +127,14 @@ Mesh::Mesh( const MEDCoupling::MEDCouplingIMesh* mesh )
                        originPtr+_spaceDim,
                        dxyzPtr,
                        dxyzPtr+_spaceDim);
+    _meshNotDeleted=true;
+    _isStructured=true;
+
        delete [] originPtr;
        delete [] dxyzPtr;
        delete [] nodeStrctPtr;
        delete [] Box0 ;
+
        setMesh();
 }
 
@@ -127,10 +142,10 @@ Mesh::Mesh( const MEDCoupling::MEDCouplingUMesh* mesh )
 {
        _spaceDim=mesh->getSpaceDimension();
        _meshDim=mesh->getMeshDimension();
-    _isStructured=false;
        double* Box0=new double[2*_spaceDim];
        mesh->getBoundingBox(Box0);
     _name=mesh->getName();
+    _epsilon=1e-6;
     _indexFacePeriodicSet=false;
     
        _xMin=Box0[0];
@@ -147,10 +162,21 @@ Mesh::Mesh( const MEDCoupling::MEDCouplingUMesh* mesh )
        }
 
        _mesh=mesh->deepCopy();
+       _mesh=mesh->buildUnstructured();
+    _meshNotDeleted=true;
+    _isStructured=false;
        delete [] Box0 ;
+
        setMesh();
 }
 
+//----------------------------------------------------------------------
+Mesh::Mesh( const std::string filename, int meshLevel )
+//----------------------------------------------------------------------
+{
+       readMeshMed(filename, meshLevel);
+}
+
 //----------------------------------------------------------------------
 Mesh::Mesh( const Mesh& mesh )
 //----------------------------------------------------------------------
@@ -158,6 +184,7 @@ Mesh::Mesh( const Mesh& mesh )
        _spaceDim = mesh.getSpaceDimension() ;
        _meshDim = mesh.getMeshDimension() ;
     _name=mesh.getName();
+    _epsilon=mesh.getComparisonEpsilon();
     _xMax=mesh.getXMax();
     _yMin=mesh.getYMin();
     _yMax=mesh.getYMax();
@@ -208,13 +235,7 @@ Mesh::Mesh( const Mesh& mesh )
     
        MCAuto<MEDCouplingMesh> m1=mesh.getMEDCouplingMesh()->deepCopy();
        _mesh=m1;
-}
-
-//----------------------------------------------------------------------
-Mesh::Mesh( const std::string filename, int meshLevel )
-//----------------------------------------------------------------------
-{
-       readMeshMed(filename, meshLevel);
+    _meshNotDeleted=mesh.meshNotDeleted();
 }
 
 //----------------------------------------------------------------------
@@ -229,11 +250,13 @@ Mesh::readMeshMed( const std::string filename, const int meshLevel)
        _meshDim=_mesh->getMeshDimension();
        _spaceDim=_mesh->getSpaceDimension();
     _name=_mesh->getName();
+    _epsilon=1e-6;
     _indexFacePeriodicSet=false;
     MEDCoupling::MEDCouplingIMesh* structuredMesh = dynamic_cast<MEDCoupling::MEDCouplingIMesh*> (_mesh.retn());
     if(structuredMesh)
     {
         _isStructured=true;
+               _meshNotDeleted=false;
         _dxyz=structuredMesh->getDXYZ();
         _nxyz=structuredMesh->getCellGridStructure();
         double* Box0=new double[2*_spaceDim];
@@ -253,73 +276,218 @@ Mesh::readMeshMed( const std::string filename, const int meshLevel)
         }
     }
     else
+    {
         _isStructured=false;
+               _meshNotDeleted=true;
+    }
     
        MEDCouplingUMesh*  mu = setMesh();
        setGroups(m, mu);
 
        cout<<endl<< "Loaded file "<< filename<<endl;
-    cout<<"Mesh name= "<<m->getName()<<", mesh dim="<< _meshDim<< ", space dim="<< _spaceDim<< ", nb cells= "<<getNumberOfCells()<< ", nb nodes= "<<getNumberOfNodes()<<endl;
+    cout<<"Mesh name = "<<m->getName()<<", mesh dim = "<< _meshDim<< ", space dim = "<< _spaceDim<< ", nb cells= "<<getNumberOfCells()<< ", nb nodes= "<<getNumberOfNodes()<<endl;
 
        m->decrRef();
-       mu->decrRef();
 }
 
 void
-Mesh::setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName)
+Mesh::setGroupAtFaceByCoords(double x, double y, double z, double eps, std::string groupName, bool isBoundaryGroup)
 {
-       int nbFace=getNumberOfFaces();
-       bool flag=false;
-       for (int iface=0;iface<nbFace;iface++)
-       {
-               double FX=_faces[iface].x();
-               double FY=_faces[iface].y();
-               double FZ=_faces[iface].z();
-               if (abs(FX-x)<eps && abs(FY-y)<eps && abs(FZ-z)<eps)
+       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++)
                {
-                       _faces[iface].setGroupName(groupName);
-                       std::vector< int > nodesID= _faces[iface].getNodesId();
-                       int nbNodes = _faces[iface].getNumberOfNodes();
-                       for(int inode=0 ; inode<nbNodes ; inode++)
-                               _nodes[nodesID[inode]].setGroupName(groupName);
+                       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);
+                       }
+               }
 
-                       flag=true;
+       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 (flag)
+}
+
+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[_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 (nodeIds.size()>0)
     {
-               _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
-               _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
-        //To do : update _faceGroups and _nodeGroups
+               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)
+Mesh::setGroupAtPlan(double value, int direction, double eps, std::string groupName, bool isBoundaryGroup)
 {
-       int nbFace=getNumberOfFaces();
-       bool flag=false;
-       for (int iface=0;iface<nbFace;iface++)
-       {
-               double cord=_faces[iface].getBarryCenter()[direction];
-               if (abs(cord-value)<eps)
+       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++)
                {
-                       _faces[iface].setGroupName(groupName);
-                       std::vector< int > nodesID= _faces[iface].getNodesId();
-                       int nbNodes = _faces[iface].getNumberOfNodes();
-                       for(int inode=0 ; inode<nbNodes ; inode++)
-                {
-                               _nodes[nodesID[inode]].setGroupName(groupName);
-                }
+                       cord=_faces[iface].getBarryCenter()[direction];
+                       if (abs(cord-value)<eps)
+                       {
+                               faceIds.insert(faceIds.end(),iface);
+                               _faces[iface].setGroupName(groupName);
+                       }
+               }
 
-                       flag=true;
+       /* 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);
+                       }
+               }
+
+       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 (flag)
+       if (nodeIds.size()>0)
     {
-               _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
-               _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
-        //To do : update _faceGroups, _nodeGroups
+               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;
+               }
     }
 }
 
@@ -351,8 +519,6 @@ Mesh::setPeriodicFaces(bool check_groups, bool use_central_inversion)
     if(_indexFacePeriodicSet)
         return;
         
-    double eps=1.E-10;
-
     for (int indexFace=0;indexFace<_boundaryFaceIds.size() ; indexFace++)
     {
         Face my_face=_faces[_boundaryFaceIds[indexFace]];
@@ -376,13 +542,13 @@ Mesh::setPeriodicFaces(bool check_groups, bool use_central_inversion)
                 Face face_i=_faces[_boundaryFaceIds[iface]];
                 double xi=face_i.x();
                 double yi=face_i.y();
-                if (   (abs(y-yi)<eps || abs(x-xi)<eps )// Case of a square geometry
+                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)<eps ) // Case of a central inversion
-                    && fabs(my_face.getMeasure() - face_i.getMeasure())<eps
-                    && fabs(my_face.getXN()      + face_i.getXN())<eps
-                    && fabs(my_face.getYN()      + face_i.getYN())<eps
-                    && fabs(my_face.getZN()      + face_i.getZN())<eps )
+                    && ( !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;
@@ -401,13 +567,13 @@ Mesh::setPeriodicFaces(bool check_groups, bool use_central_inversion)
                 double xi=face_i.x();
                 double yi=face_i.y();
                 double zi=face_i.z();
-                if ( ((abs(y-yi)<eps && abs(x-xi)<eps) || (abs(x-xi)<eps && abs(z-zi)<eps) || (abs(y-yi)<eps && abs(z-zi)<eps))// Case of a cube geometry
+                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)<eps )// Case of a central inversion
-                    && fabs(my_face.getMeasure() - face_i.getMeasure())<eps
-                    && fabs(my_face.getXN()      + face_i.getXN())<eps
-                    && fabs(my_face.getYN()      + face_i.getYN())<eps
-                    && fabs(my_face.getZN()      + face_i.getZN())<eps )
+                    && ( !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;
@@ -542,7 +708,10 @@ Mesh::getElementTypesNames() const
 void
 Mesh::setGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh*  mu)
 {
-       //Searching for face groups
+       int nbCellsSubMesh, nbNodesSubMesh;
+       bool foundFace, foundNode;
+       
+       /* Searching for face groups */
        vector<string> faceGroups=medmesh->getGroupsNames() ;
 
        for (unsigned int i=0;i<faceGroups.size();i++ )
@@ -554,15 +723,18 @@ Mesh::setGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh*  mu)
                vector<mcIdType>::iterator it = find(nonEmptyGrp.begin(), nonEmptyGrp.end(), -1);
                if (it != nonEmptyGrp.end())
                {
-                       cout<<"Boundary face group named "<< groupName << " found"<<endl;
                        MEDCouplingUMesh *m=medmesh->getGroup(-1,groupName.c_str());
-            mu->unPolyze();
-                       _faceGroups.insert(_faceGroups.begin(),m);
-                       _faceGroupNames.insert(_faceGroupNames.begin(),groupName);
+            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();
-
-                       int nbCellsSubMesh=m->getNumberOfCells();
+                       /* Face identification */
                        for (int ic(0), k(0); ic<nbCellsSubMesh; ic++, k+=_spaceDim)
                        {
                                vector<double> coorBaryXyz(3,0);
@@ -570,41 +742,27 @@ Mesh::setGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh*  mu)
                                        coorBaryXyz[d] = coorBary[k+d];
                                Point p1(coorBaryXyz[0],coorBaryXyz[1],coorBaryXyz[2]) ;
 
-                               int flag=0;
-                /* double min=INFINITY;
-                Point p3;
-                int closestFaceNb;*/
+                               foundFace=false;
                                for (int iface=0;iface<_numberOfFaces;iface++ )
                                {
                                        Point p2=_faces[iface].getBarryCenter();
-                    /*if(groupName=="Wall" and p1.distance(p2)<min)
-                    {
-                        min=p1.distance(p2);
-                        p3=p2;
-                        closestFaceNb=iface;
-                    }*/
-                                       if(p1.distance(p2)<1.E-10)
+                                       if(p1.distance(p2)<_epsilon)
                                        {
                                                _faces[iface].setGroupName(groupName);
-                                               flag=1;
+                                               _faceGroupsIds[_faceGroupsIds.size()-1][ic]=iface;
+                                               foundFace=true;
                                                break;
                                        }
                                }
-                /* if(groupName=="Wall" and min>1.E-10)
-                    {
-                        cout.precision(15);
-                        cout<<"Subcell number "<< ic <<" Min distance to Wall = "<<min <<" p= "<<p1[0]<<" , "<<p1[1]<<" , "<<p1[2]<<endl;
-                        cout<<" attained at face "<< closestFaceNb << " p_min= "<<p3[0]<<" , "<<p3[1]<<" , "<<p3[2]<<endl;
-                    }*/
-                               if (flag==0)
-                                       throw CdmathException("No face belonging to group " + groupName + " found");
+                               if (not foundFace)
+                                       throw CdmathException("No face found for group " + groupName );
                        }
                        baryCell->decrRef();
                        //m->decrRef();
                }
        }
 
-       //Searching for node groups
+       /* Searching for node groups */
        vector<string> nodeGroups=medmesh->getGroupsOnSpecifiedLev(1) ;
 
        for (unsigned int i=0;i<nodeGroups.size();i++ )
@@ -615,16 +773,16 @@ Mesh::setGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh*  mu)
 
                if(nodeids!=NULL)
                {
-                       cout<<"Boundary node group named "<< groupName << " found"<<endl;
+                       _nodeGroups.insert(_nodeGroups.end(),nodeGroup);
+                       _nodeGroupNames.insert(_nodeGroupNames.end(),groupName);
 
-                       _nodeGroups.insert(_nodeGroups.begin(),nodeGroup);
-                       _nodeGroupNames.insert(_nodeGroupNames.begin(),groupName);
-
-                       int nbNodesSubMesh=nodeGroup->getNumberOfTuples();//nodeGroup->getNbOfElems();
+                       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);
@@ -632,19 +790,20 @@ Mesh::setGroups( const MEDFileUMesh* medmesh, MEDCouplingUMesh*  mu)
                                        coorP[d] = cood[nodeids[ic]*_spaceDim+d];
                                Point p1(coorP[0],coorP[1],coorP[2]) ;
 
-                               int flag=0;
+                               foundNode=false;
                                for (int inode=0;inode<_numberOfNodes;inode++ )
                                {
                                        Point p2=_nodes[inode].getPoint();
-                                       if(p1.distance(p2)<1.E-10)
+                                       if(p1.distance(p2)<_epsilon)
                                        {
                                                _nodes[inode].setGroupName(groupName);
-                                               flag=1;
+                                               _nodeGroupsIds[_nodeGroupsIds.size()-1][ic]=inode;
+                                               foundNode=true;
                                                break;
                                        }
                                }
-                               if (flag==0)
-                                       throw CdmathException("No node belonging to group " + groupName + " found");
+                               if (not foundNode)
+                                       throw CdmathException("No node found for group " + groupName );
                        }
                }
        }
@@ -655,26 +814,34 @@ MEDCouplingUMesh*
 Mesh::setMesh( void )
 //----------------------------------------------------------------------
 {
+       /* This is the main function translating medcouplingumesh info into Mesh class to be used when designing numerical methods
+        * We need the level 0 mesh to extract the cell-node connectvity
+        * We need the level -1 mesh to extract the cell-face and face-node connectivities (use o build descending connectivity)
+        * Be careful : the nodes in the medcoupling mesh are not necessarily all conected to a cell/face. 
+        * Mesh class discard isolated nodes, hence the number of nodes in Mesh class can be lower than the number of nodes in medcouplingumesh.
+        */
+        
        DataArrayIdType *desc  = DataArrayIdType::New();
        DataArrayIdType *descI = DataArrayIdType::New();
        DataArrayIdType *revDesc  = DataArrayIdType::New();
        DataArrayIdType *revDescI = DataArrayIdType::New();
        MEDCouplingUMesh* mu = _mesh->buildUnstructured();
-
+       MEDCouplingUMesh* mu2;//mesh of dimension N-1 containing the cell interfaces->cell/face connectivity
+       
        mu->unPolyze();
-       /* Save the boundary mesh for later use*/
-       _boundaryMesh = mu->computeSkin();
+    mu->sortCellsInMEDFileFrmt( );
        
-       MEDCouplingUMesh* mu2=mu->buildDescendingConnectivity2(desc,descI,revDesc,revDescI);//mesh of dimension N-1 containing the cell interfaces
-    
-    const mcIdType *tmp = desc->getConstPointer();
+       if(_meshDim<2)
+               mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
+       else
+               mu2=mu->buildDescendingConnectivity2(desc,descI,revDesc,revDescI);
+       
+    const mcIdType *tmp = desc->getConstPointer();//Lists the faces surrounding each cell
     const mcIdType *tmpI=descI->getConstPointer();
 
-       const mcIdType *tmpA =revDesc->getConstPointer();
+       const mcIdType *tmpA =revDesc->getConstPointer();//Lists the cells surrounding each face
        const mcIdType *tmpAI=revDescI->getConstPointer();
 
-    //const int *work=tmp+tmpI[id];//corresponds to buildDescendingConnectivity
-
        //Test du type d'éléments contenus dans le maillage afin d'éviter les éléments contenant des points de gauss
        _eltsTypes=mu->getAllGeoTypesSorted();
        for(int i=0; i<_eltsTypes.size();i++)
@@ -694,37 +861,37 @@ Mesh::setMesh( void )
        }
 
        DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
-       const double *coorBary=baryCell->getConstPointer();
+       const double *coorBary=baryCell->getConstPointer();//Used for cell center coordinates
 
        MEDCouplingFieldDouble* fields=mu->getMeasureField(true);
        DataArrayDouble *surface = fields->getArray();
-       const double *surf=surface->getConstPointer();
+       const double *surf=surface->getConstPointer();//Used for cell lenght/surface/volume
 
        DataArrayDouble *coo = mu->getCoords() ;
-       const double    *cood=coo->getConstPointer();
+       const double    *cood=coo->getConstPointer();//Used for nodes coordinates
 
        DataArrayIdType *revNode =DataArrayIdType::New();
        DataArrayIdType *revNodeI=DataArrayIdType::New();
        mu->getReverseNodalConnectivity(revNode,revNodeI) ;
-       const mcIdType *tmpN =revNode->getConstPointer();
+       const mcIdType *tmpN =revNode->getConstPointer();//Used to know which cells surround a given node
        const mcIdType *tmpNI=revNodeI->getConstPointer();
 
        DataArrayIdType *revCell =DataArrayIdType::New();
        DataArrayIdType *revCellI=DataArrayIdType::New();
-       mu2->getReverseNodalConnectivity(revCell,revCellI) ;
-       const mcIdType *tmpC =revCell->getConstPointer();
+       mu2->getReverseNodalConnectivity(revCell,revCellI);
+       const mcIdType *tmpC =revCell->getConstPointer();//Used to know which faces surround a given node
        const mcIdType *tmpCI=revCellI->getConstPointer();
 
        const DataArrayIdType *nodal  = mu2->getNodalConnectivity() ;
        const DataArrayIdType *nodalI = mu2->getNodalConnectivityIndex() ;
-       const mcIdType *tmpNE =nodal->getConstPointer();
+       const mcIdType *tmpNE =nodal->getConstPointer();//Used to know which nodes surround a given face
        const mcIdType *tmpNEI=nodalI->getConstPointer();
 
        _numberOfCells = mu->getNumberOfCells() ;
        _cells      = new Cell[_numberOfCells] ;
 
-       _numberOfNodes = mu->getNumberOfNodes() ;
-       _nodes      = new Node[_numberOfNodes] ;
+       _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
 
        _numberOfFaces = mu2->getNumberOfCells();
        _faces       = new Face[_numberOfFaces] ;
@@ -756,86 +923,134 @@ Mesh::setMesh( void )
     }    
 
        // _cells, _nodes and _faces initialization:
-       if (_spaceDim == 1)
+       if (_meshDim == 1)
        {
+               double xn, yn=0., zn=0.;//Components of the normal vector at a cell interface
+               double norm;
                for( int id=0;id<_numberOfCells;id++ )
                {
-                       const mcIdType *work=tmp+tmpI[id];
-            int nbFaces=tmpI[id+1]-tmpI[id];
-                       
-                       int nbVertices=mu->getNumberOfNodesInCell(id) ;
-            
-                       Cell ci( nbVertices, nbFaces, surf[id], Point(coorBary[id], 0.0, 0.0) ) ;
-
+                       Point p(0.0,0.0,0.0) ;
+                       for(int idim=0; idim<_spaceDim; idim++)
+                               p[idim]=coorBary[id*_spaceDim+idim];
+       
+                       mcIdType nbVertices=mu->getNumberOfNodesInCell(id) ;//should be equal to 2
+                       assert( nbVertices==2);
                        std::vector<mcIdType> nodeIdsOfCell ;
                        mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
-                       for( int el=0;el<nbVertices;el++ )
-            {
-                               ci.addNodeId(el,nodeIdsOfCell[el]) ;
-                ci.addFaceId(el,nodeIdsOfCell[el]) ;
-            }
-            
-            //This assumes that in a 1D mesh the cell numbers form a consecutive sequence 
-            //going either from left to right (xn=-1) or right to left (xn=1)
-                       double xn = (cood[nodeIdsOfCell[nbVertices-1]] - cood[nodeIdsOfCell[0]] > 0.0) ? -1.0 : 1.0;
-
-                       for( int el=0;el<nbFaces;el++ )
+       
+               mcIdType nbFaces=tmpI[id+1]-tmpI[id];//should be equal to 2
+                       assert( nbFaces==2);
+               const mcIdType *work=tmp+tmpI[id];
+       
+                       /* compute the normal to the face */
+                   xn = cood[nodeIdsOfCell[0]*_spaceDim  ] - cood[nodeIdsOfCell[nbFaces-1]*_spaceDim  ];
+               if(_spaceDim>1)        
+                               yn = cood[nodeIdsOfCell[0]*_spaceDim+1] - cood[nodeIdsOfCell[nbFaces-1]*_spaceDim+1];
+               if(_spaceDim>2)        
+                               zn = cood[nodeIdsOfCell[0]*_spaceDim+2] - cood[nodeIdsOfCell[nbFaces-1]*_spaceDim+2];
+                       norm = sqrt(xn*xn+yn*yn+zn*zn);
+                       if(norm<_epsilon)
+                               throw CdmathException("!!! Mesh::setMesh Normal vector has norm 0 !!!");
+                       else
                        {
-                               ci.addNormalVector(el,xn,0.0,0.0) ;
-                               int indexFace=abs(work[el])-1;
-                ci.addFaceId(el,indexFace) ;
-                               xn = - xn;
+                               xn /= norm;
+                               yn /= norm;
+                               zn /= norm;
+                       }
+               
+                       Cell ci( nbVertices, nbFaces, surf[id], p ) ;//nbCells=nbFaces=2
+               for( int el=0;el<nbFaces;el++ )
+                       {
+                               ci.addNodeId(el,nodeIdsOfCell[el]) ;//global node number
+                               ci.addNormalVector(el,xn,yn,zn) ;
+                               ci.addFaceId(el,work[el]) ;
+                               xn = - xn; yn=-yn; zn=-zn;
                        }
                        _cells[id] = ci ;
                }
-
-               for( int id(0), k(0); id<_numberOfNodes; id++, k+=_spaceDim)
-               {
-                       Point p(cood[k], 0.0, 0.0) ;
-                       const mcIdType *workc=tmpN+tmpNI[id];
-                       mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
-
-                       const mcIdType *workf=tmpC+tmpCI[id];
-                       mcIdType nbFaces=tmpCI[id+1]-tmpCI[id];
-            const mcIdType *workn=tmpA+tmpAI[id];
-            mcIdType nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
-                       Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
-
-                       for( int el=0;el<nbCells;el++ )
-                               vi.addCellId(el,workc[el]) ;
-                       for( int el=0;el<nbFaces;el++ )
-                               vi.addFaceId(el,workf[el]) ;
-                       for( int el=0;el<nbNeighbourNodes;el++ )
-                               vi.addNeighbourNodeId(el,workn[el]) ;
-                       _nodes[id] = vi ;
-               }
-        _boundaryNodeIds.push_back(0);
-        _boundaryNodeIds.push_back(_numberOfNodes-1);
-
-               for(int id(0), k(0); id<_numberOfFaces; id++, k+=_spaceDim)
+       
+               for( int id(0); id<_numberOfFaces; id++ )
                {
-                       Point p(cood[k], 0.0, 0.0) ;
-                       const mcIdType *workc=tmpA+tmpAI[id];
-                       mcIdType nbCells=tmpAI[id+1]-tmpAI[id];
-
                        const mcIdType *workv=tmpNE+tmpNEI[id]+1;
-                       Face fi( 1, nbCells, 1.0, p, 1.0, 0.0, 0.0) ;
-                       fi.addNodeId(0,workv[0]) ;
-
-                       for(int idCell=0; idCell<nbCells; idCell++)
-                               fi.addCellId(idCell,workc[idCell]) ;
-
+                       mcIdType nbNodes= tmpNEI[id+1]-tmpNEI[id]-1;//Normally equal to 1.
+                       assert( nbNodes==1);
+       
+                       std::vector<double> coo(0) ;
+                       mu2->getCoordinatesOfNode(workv[0],coo);
+                       Point p(0,0.0,0.0) ;
+                       for(int idim=0; idim<_spaceDim; idim++)
+                               p[idim]=coo[idim];
+       
+                   const mcIdType *workc=tmpA+tmpAI[id];
+                   mcIdType nbCells=tmpAI[id+1]-tmpAI[id];
+                       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
+                               fi.addNodeId(node_id,workv[node_id]) ;//global node number
+       
+                       fi.addCellId(0,workc[0]) ;
+                       for(int cell_id=1; cell_id<nbCells;cell_id++)
+                       {
+                               int cell_idx=0;
+                               if (workc[cell_id]!=workc[cell_id-1])//For some meshes (bad ones) the same cell can appear several times
+                                       {
+                                       fi.addCellId(cell_idx+1,workc[cell_id]) ;
+                                       cell_idx++;
+                                       }                
+                       }
+                       if(nbCells==1)
+                               _boundaryFaceIds.push_back(id);
                        _faces[id] = fi ;
                }
-        _boundaryFaceIds.push_back(0);
-        _boundaryFaceIds.push_back(_numberOfFaces-1);
+       
+               int correctNbNodes=0;
+               for( int id=0;id<_numberOfNodes;id++ )
+               {
+                       const mcIdType *workc=tmpN+tmpNI[id];
+                       mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
+                       
+                       if( nbCells>0)//To make sure this is not an isolated node
+                       {
+                               correctNbNodes++;
+                               std::vector<double> coo(0) ;
+                               mu->getCoordinatesOfNode(id,coo);
+                               Point p(0,0.0,0.0) ;
+                               for(int idim=0; idim<_spaceDim; idim++)
+                                       p[idim]=coo[idim];
+               
+                               const mcIdType *workf=tmpC+tmpCI[id];
+                               mcIdType nbFaces=tmpCI[id+1]-tmpCI[id];
+                               assert( nbFaces==1);
+               
+                           const mcIdType *workn=tmpN+tmpNI[id];
+                           mcIdType nbNeighbourNodes=tmpNI[id+1]-tmpNI[id];
+                       
+                               Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
+                       for( int el=0;el<nbCells;el++ )
+                                       vi.addCellId(el,workc[el]) ;
+                       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()) ;
+                               if(vi.isBorder())
+                                       _boundaryNodeIds.push_back(id);
+                               _nodes[id] = vi ;
+                       }
+               }
+               if( _numberOfNodes!=correctNbNodes)
+                       cout<<"Found isolated nodes : correctNbNodes= "<<correctNbNodes<<", _numberOfNodes= "<<_numberOfNodes<<endl;
        }
-       else if(_spaceDim==2  || _spaceDim==3)
+       else if(_meshDim==2  || _meshDim==3)
        {
-               DataArrayDouble *barySeg = mu2->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;
+               DataArrayDouble *barySeg = mu2->computeIsoBarycenterOfNodesPerCell();//computeCellCenterOfMass() ;//Used as face center
                const double *coorBarySeg=barySeg->getConstPointer();
 
-               MEDCouplingFieldDouble* fieldn;
+               MEDCouplingFieldDouble* fieldl=mu2->getMeasureField(true);
+               DataArrayDouble *longueur = fieldl->getArray();
+               const double *lon=longueur->getConstPointer();//The lenght/area of each face
+
+               MEDCouplingFieldDouble* fieldn;//The normal to each face
                DataArrayDouble *normal;
                const double *tmpNormal;
 
@@ -884,7 +1099,7 @@ Mesh::setMesh( void )
                                        ci.addFaceId(el,faceIndex) ;
                                }
                        else//build normals associated to the couple (cell id, face el)
-                       {
+                       {//Case _meshDim=1 should be moved upper since we are in the 2D/3D branch
                                if(_meshDim==1)//we know in this case there are only two faces around the cell id, each face is composed of a single node
                                {//work[0]= first face global number, work[1]= second face global number
                     mcIdType indexFace0=abs(work[0])-1;//=work[0] since Fortran type numbering was used, and negative sign means anticlockwise numbering
@@ -948,10 +1163,6 @@ Mesh::setMesh( void )
                        _cells[id] = ci ;
                }
 
-               MEDCouplingFieldDouble* fieldl=mu2->getMeasureField(true);
-               DataArrayDouble *longueur = fieldl->getArray();
-               const double *lon=longueur->getConstPointer();
-
                if(_spaceDim!=_meshDim)
                {
                        /* Since spaceDim!=meshDim, don't build normal to faces */
@@ -963,7 +1174,7 @@ Mesh::setMesh( void )
                /*Building mesh faces */
                for(int id(0), k(0); id<_numberOfFaces; id++, k+=_spaceDim)
                {
-                       vector<double> coorBarySegXyz(3,0);
+                       vector<double> coorBarySegXyz(3);
                        for (int d=0; d<_spaceDim; d++)
                                coorBarySegXyz[d] = coorBarySeg[k+d];
                        Point p(coorBarySegXyz[0],coorBarySegXyz[1],coorBarySegXyz[2]) ;
@@ -1012,51 +1223,55 @@ Mesh::setMesh( void )
                }
 
                /*Building mesh nodes, should be done after building mesh faces in order to detect boundary nodes*/
+               int correctNbNodes=0;
                for(int id(0), k(0); id<_numberOfNodes; id++, k+=_spaceDim)
                {
-                       vector<double> coorP(3,0);
-                       for (int d=0; d<_spaceDim; d++)
-                               coorP[d] = cood[k+d];
-                       Point p(coorP[0],coorP[1],coorP[2]) ;
-
                        const mcIdType *workc=tmpN+tmpNI[id];
                        mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
-                       const mcIdType *workf=tmpC+tmpCI[id];
-                       mcIdType nbFaces=tmpCI[id+1]-tmpCI[id];
-                       const mcIdType *workn;
-                       mcIdType nbNeighbourNodes;
-            if (_meshDim == 1)
-            {
-                workn=tmpA+tmpAI[id];
-                nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
-            }
-            else if (_meshDim == 2)
-            {
-                workn=tmpC+tmpCI[id];
-                nbNeighbourNodes=tmpCI[id+1]-tmpCI[id];
-            }
-            else//_meshDim == 3
-            {
-                workn=tmpN2+tmpNI2[id];
-                nbNeighbourNodes=tmpNI2[id+1]-tmpNI2[id];
-            }    
-                       Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
-
-                       for( int el=0;el<nbCells;el++ )
-                               vi.addCellId(el,workc[el]) ;
-                       for( int el=0;el<nbNeighbourNodes;el++ )
-                               vi.addNeighbourNodeId(el,workn[el]) ;
-            //Detection of border nodes    
-            bool isBorder=false;
-                       for( int el=0;el<nbFaces;el++ )
-            {
-                               vi.addFaceId(el,workf[el],_faces[workf[el]].isBorder()) ;
-                isBorder= isBorder || _faces[workf[el]].isBorder();
-            }
-            if(isBorder)
-                _boundaryNodeIds.push_back(id);
-                       _nodes[id] = vi ;
+
+                       if( nbCells>0)//To make sure this is not an isolated node
+                       {
+                               correctNbNodes++;
+                               vector<double> coorP(3);
+                               for (int d=0; d<_spaceDim; d++)
+                                       coorP[d] = cood[k+d];
+                               Point p(coorP[0],coorP[1],coorP[2]) ;
+               
+                               const mcIdType *workf=tmpC+tmpCI[id];
+                               mcIdType nbFaces=tmpCI[id+1]-tmpCI[id];
+                               const mcIdType *workn;
+                               mcIdType nbNeighbourNodes;
+                               if (_meshDim == 1)
+                               {
+                                       workn=tmpA+tmpAI[id];
+                                       nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
+                               }
+                               else if (_meshDim == 2)
+                               {
+                                       workn=tmpC+tmpCI[id];
+                                       nbNeighbourNodes=tmpCI[id+1]-tmpCI[id];
+                               }
+                               else//_meshDim == 3
+                               {
+                                       workn=tmpN2+tmpNI2[id];
+                                       nbNeighbourNodes=tmpNI2[id+1]-tmpNI2[id];
+                               }    
+                               Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
+               
+                               for( int el=0;el<nbCells;el++ )
+                                       vi.addCellId(el,workc[el]) ;
+                               for( int el=0;el<nbNeighbourNodes;el++ )
+                                       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()) ;
+                               if(vi.isBorder())
+                                       _boundaryNodeIds.push_back(id);
+                               _nodes[id] = vi ;
+                       }
                }
+               if( _numberOfNodes!=correctNbNodes)
+                       cout<<"Found isolated nodes : correctNbNodes= "<<correctNbNodes<<", _numberOfNodes= "<<_numberOfNodes<<endl;
 
                if(_spaceDim==_meshDim)
                        fieldn->decrRef();
@@ -1085,11 +1300,21 @@ Mesh::setMesh( void )
             _zMax=Box0[5];
         }
     }
-
+       
     //Set boundary groups
-       _faceGroupNames.push_back("Boundary");
+    _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();
@@ -1112,185 +1337,10 @@ Mesh::setMesh( void )
         revDescI2->decrRef();
         mu3->decrRef();
     }
-       
+       
     return mu;
 }
 
-//----------------------------------------------------------------------
-Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName )
-//----------------------------------------------------------------------
-{
-       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;
-    _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);
-       delete [] originPtr;
-       delete [] dxyzPtr;
-       delete [] nodeStrctPtr;
-
-       DataArrayIdType *desc=DataArrayIdType::New();
-       DataArrayIdType *descI=DataArrayIdType::New();
-       DataArrayIdType *revDesc=DataArrayIdType::New();
-       DataArrayIdType *revDescI=DataArrayIdType::New();
-       MEDCouplingUMesh* mu=_mesh->buildUnstructured();
-       MEDCouplingUMesh *mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
-
-       const mcIdType *tmp=desc->getConstPointer();
-       const mcIdType *tmpI=descI->getConstPointer();
-
-       const mcIdType *tmpA =revDesc->getConstPointer();
-       const mcIdType *tmpAI=revDescI->getConstPointer();
-
-       _eltsTypes=mu->getAllGeoTypesSorted();
-
-       DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
-       const double *coorBary=baryCell->getConstPointer();
-
-       _numberOfCells = _mesh->getNumberOfCells() ;
-       _cells    = new Cell[_numberOfCells] ;
-
-       _numberOfNodes = mu->getNumberOfNodes() ;
-       _nodes    = new Node[_numberOfNodes] ;
-
-       _numberOfFaces = _numberOfNodes;
-       _faces    = new Face[_numberOfFaces] ;
-
-    _numberOfEdges = _numberOfCells;
-    
-       MEDCouplingFieldDouble* fieldl=mu->getMeasureField(true);
-       DataArrayDouble *longueur = fieldl->getArray();
-       const double *lon=longueur->getConstPointer();
-
-       DataArrayDouble *coo = mu->getCoords() ;
-       const double *cood=coo->getConstPointer();
-
-       int comp=0;
-       for( int id=0;id<_numberOfCells;id++ )
-       {
-               int nbVertices=mu->getNumberOfNodesInCell(id) ;
-               Point p(coorBary[id],0.0,0.0) ;
-               Cell ci( nbVertices, nbVertices, lon[id], p ) ;
-
-               std::vector<mcIdType> nodeIdsOfCell ;
-               mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
-               for( int el=0;el<nbVertices;el++ )
-               {
-                       ci.addNodeId(el,nodeIdsOfCell[el]) ;
-                       ci.addFaceId(el,nodeIdsOfCell[el]) ;
-               }
-
-        double xn = (cood[nodeIdsOfCell[nbVertices-1]] - cood[nodeIdsOfCell[0]] > 0.0) ? -1.0 : 1.0;
-
-        mcIdType nbFaces=tmpI[id+1]-tmpI[id];
-        const mcIdType *work=tmp+tmpI[id];
-               
-        for( int el=0;el<nbFaces;el++ )
-               {
-                       ci.addNormalVector(el,xn,0.0,0.0) ;
-                       ci.addFaceId(el,work[el]) ;
-                       xn = - xn;
-               }
-
-               _cells[id] = ci ;
-
-               comp=comp+2;
-       }
-
-    //Suppress the following since tmpN=tmpA
-       DataArrayIdType *revNode=DataArrayIdType::New();
-       DataArrayIdType *revNodeI=DataArrayIdType::New();
-       mu->getReverseNodalConnectivity(revNode,revNodeI) ;
-       const mcIdType *tmpN=revNode->getConstPointer();
-       const mcIdType *tmpNI=revNodeI->getConstPointer();
-
-       for( int id=0;id<_numberOfNodes;id++ )
-       {
-               std::vector<double> coo ;
-               mu->getCoordinatesOfNode(id,coo);
-               Point p(coo[0],0.0,0.0) ;
-               const mcIdType *workc=tmpN+tmpNI[id];
-               mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
-               mcIdType nbFaces=1;
-        const mcIdType *workn=tmpA+tmpAI[id];
-        mcIdType nbNeighbourNodes=tmpAI[id+1]-tmpAI[id];
-        
-               Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
-        for( int el=0;el<nbCells;el++ )
-                       vi.addCellId(el,workc[el]) ;
-               for( int el=0;el<nbFaces;el++ )
-                       vi.addFaceId(el,id) ;
-        for( int el=0;el<nbNeighbourNodes;el++ )
-                       vi.addNeighbourNodeId(el,workn[el]) ;
-               _nodes[id] = vi ;
-
-
-               int nbVertices=1;
-               Face fi( nbVertices, nbCells, 1.0, p, 1., 0., 0. ) ;
-        for( int el=0;el<nbVertices;el++ )
-                       fi.addNodeId(el,id) ;
-
-               for( int el=0;el<nbCells;el++ )
-                       fi.addCellId(el,workc[el]) ;
-               _faces[id] = fi ;
-       }
-
-    //Set boundary groups
-       _faceGroupNames.push_back("Boundary");
-    _nodeGroupNames.push_back("Boundary");
-    _boundaryFaceIds.push_back(0);
-    _boundaryFaceIds.push_back(_numberOfFaces-1);
-    _boundaryNodeIds.push_back(0);
-    _boundaryNodeIds.push_back(_numberOfNodes-1);
-
-       fieldl->decrRef();
-       baryCell->decrRef();
-       desc->decrRef();
-       descI->decrRef();
-       revDesc->decrRef();
-       revDescI->decrRef();
-       revNode->decrRef();
-       revNodeI->decrRef();
-       mu2->decrRef();
-       mu->decrRef();
-}
-
 //----------------------------------------------------------------------
 Mesh::Mesh( std::vector<double> points, std::string meshName )
 //----------------------------------------------------------------------
@@ -1311,6 +1361,7 @@ Mesh::Mesh( std::vector<double> points, std::string meshName )
        _spaceDim = 1 ;
        _meshDim  = 1 ;
     _name=meshName;
+    _epsilon=1e-6;
        _xMin=points[0];
        _xMax=points[nx-1];
        _yMin=0.;
@@ -1343,110 +1394,67 @@ Mesh::Mesh( std::vector<double> points, std::string meshName )
     delete [] coords, nodal_con;
     coords_arr->decrRef();
 
-    _mesh=mesh1d;
-    
-       DataArrayIdType *desc=DataArrayIdType::New();
-       DataArrayIdType *descI=DataArrayIdType::New();
-       DataArrayIdType *revDesc=DataArrayIdType::New();
-       DataArrayIdType *revDescI=DataArrayIdType::New();
-    MEDCouplingUMesh* mu=_mesh->buildUnstructured();
-       MEDCouplingUMesh *mu2=mu->buildDescendingConnectivity(desc,descI,revDesc,revDescI);
-
-       DataArrayDouble *baryCell = mu->computeCellCenterOfMass() ;
-       const double *coorBary=baryCell->getConstPointer();
+    _mesh=mesh1d->buildUnstructured();//To enable writeMED. Because we declared the mesh as unstructured, we decide to build the unstructured data (not mandatory)
+    _meshNotDeleted=true;
 
-       _numberOfCells = _mesh->getNumberOfCells() ;
-       _cells    = new Cell[_numberOfCells] ;
+       setMesh();
+}
 
-    _numberOfEdges = _numberOfCells;
+//----------------------------------------------------------------------
+Mesh::Mesh( double xmin, double xmax, int nx, std::string meshName )
+//----------------------------------------------------------------------
+{
+       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");
 
-       _eltsTypes=mu->getAllGeoTypesSorted();
+       double dx = (xmax - xmin)/nx ;
 
-       MEDCouplingFieldDouble* fieldl=mu->getMeasureField(true);
-       DataArrayDouble *longueur = fieldl->getArray();
-       const double *lon=longueur->getConstPointer();
+       _spaceDim = 1 ;
+       _meshDim  = 1 ;
+    _name=meshName;
+    _epsilon=1e-6;
+    _isStructured = true;
+    _indexFacePeriodicSet=false;
 
-       int comp=0;
-       for( int id=0;id<_numberOfCells;id++ )
-       {
-               int nbVertices=mu->getNumberOfNodesInCell(id) ;
-               Point p(coorBary[id],0.0,0.0) ;
-               Cell ci( nbVertices, nbVertices, lon[id], p ) ;
+       _xMin=xmin;
+       _xMax=xmax;
+       _yMin=0.;
+       _yMax=0.;
+       _zMin=0.;
+       _zMax=0.;
 
-               std::vector<mcIdType> nodeIdsOfCell ;
-               mu->getNodeIdsOfCell(id,nodeIdsOfCell) ;
-               for( int el=0;el<nbVertices;el++ )
-               {
-                       ci.addNodeId(el,nodeIdsOfCell[el]) ;
-                       ci.addFaceId(el,nodeIdsOfCell[el]) ;
-               }
-               _cells[id] = ci ;
-               comp=comp+2;
-       }
+       _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];
 
-    //Suppress the following since tmpN=tmpA
-       DataArrayIdType *revNode=DataArrayIdType::New();
-       DataArrayIdType *revNodeI=DataArrayIdType::New();
-       mu->getReverseNodalConnectivity(revNode,revNodeI) ;
-       const mcIdType *tmpN=revNode->getConstPointer();
-       const mcIdType *tmpNI=revNodeI->getConstPointer();
+       originPtr[0]=xmin;
+       nodeStrctPtr[0]=nx+1;
+       dxyzPtr[0]=dx;
 
-       _numberOfNodes = mu->getNumberOfNodes() ;
-       _nodes    = new Node[_numberOfNodes] ;
-       _numberOfFaces = _numberOfNodes;
-       _faces    = new Face[_numberOfFaces] ;
+       _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
 
-       for( int id=0;id<_numberOfNodes;id++ )
-       {
-               std::vector<double> coo ;
-               mu->getCoordinatesOfNode(id,coo);
-               Point p(coo[0],0.0,0.0) ;
-               const mcIdType *workc=tmpN+tmpNI[id];
-               mcIdType nbCells=tmpNI[id+1]-tmpNI[id];
-               mcIdType nbFaces=1;
-        const mcIdType *workn=tmpN+tmpNI[id];
-        mcIdType nbNeighbourNodes=tmpNI[id+1]-tmpNI[id];
-               Node vi( nbCells, nbFaces, nbNeighbourNodes, p ) ;
-               int nbVertices=1;
-               /* provisoire !!!!!!!!!!!!*/
-               //        Point pf(0.0,0.0,0.0) ;
-               Face fi( nbVertices, nbCells, 0.0, p, 0., 0., 0. ) ;
-
-               for( int el=0;el<nbCells;el++ )
-                       vi.addCellId(el,workc[el]) ;
-               for( int el=0;el<nbFaces;el++ )
-                       vi.addFaceId(el,id) ;
-        for( int el=0;el<nbNeighbourNodes;el++ )
-                       vi.addNeighbourNodeId(el,workn[el]) ;
-               _nodes[id] = vi ;
-
-               for( int el=0;el<nbVertices;el++ )
-                       fi.addNodeId(el,id) ;
-
-               for( int el=0;el<nbCells;el++ )
-                       fi.addCellId(el,workc[el]) ;
-               _faces[id] = fi ;
-       }
-    //Set boundary groups
-       _faceGroupNames.push_back("Boundary");
-    _nodeGroupNames.push_back("Boundary");
-    _boundaryFaceIds.push_back(0);
-    _boundaryFaceIds.push_back(_numberOfFaces-1);
-    _boundaryNodeIds.push_back(0);
-    _boundaryNodeIds.push_back(_numberOfNodes-1);
+       delete [] originPtr;
+       delete [] dxyzPtr;
+       delete [] nodeStrctPtr;
 
-       fieldl->decrRef();
-       baryCell->decrRef();
-       desc->decrRef();
-       descI->decrRef();
-       revDesc->decrRef();
-       revDescI->decrRef();
-       revNode->decrRef();
-       revNodeI->decrRef();
-       mu2->decrRef();
-       mu->decrRef();
+       setMesh();
 }
+
 //----------------------------------------------------------------------
 Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny, int split_to_triangles_policy, std::string meshName)
 //----------------------------------------------------------------------
@@ -1472,8 +1480,8 @@ Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny,
        _spaceDim = 2 ;
        _meshDim  = 2 ;
     _name=meshName;
+    _epsilon=1e-6;
     _indexFacePeriodicSet=false;
-    _isStructured = true;
        _nxyz.resize(_spaceDim);
        _nxyz[0]=nx;
        _nxyz[1]=ny;
@@ -1501,11 +1509,15 @@ Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny,
                        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;
 
     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)
         {
@@ -1536,7 +1548,7 @@ Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny,
        _spaceDim = 3;
        _meshDim  = 3;
     _name=meshName;
-    _indexFacePeriodicSet=false;
+    _epsilon=1e-6;
        _xMin=xmin;
        _xMax=xmax;
        _yMin=ymin;
@@ -1548,7 +1560,6 @@ Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny,
        double dy = (ymax - ymin)/ny ;
        double dz = (zmax - zmin)/nz ;
 
-    _isStructured = true;
        _dxyz.resize(_spaceDim);
        _dxyz[0]=dx;
        _dxyz[1]=dy;
@@ -1581,16 +1592,22 @@ Mesh::Mesh( double xmin, double xmax, int nx, double ymin, double ymax, int ny,
                        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;
 
     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 )
         {
@@ -1844,12 +1861,15 @@ Mesh::operator= ( const Mesh& mesh )
        _spaceDim = mesh.getSpaceDimension() ;
        _meshDim  = mesh.getMeshDimension() ;
     _name = mesh.getName();
+    _epsilon=mesh.getComparisonEpsilon();
        _numberOfNodes = mesh.getNumberOfNodes();
        _numberOfFaces = mesh.getNumberOfFaces();
        _numberOfCells = mesh.getNumberOfCells();
        _numberOfEdges = mesh.getNumberOfEdges();
     
     _isStructured = mesh.isStructured();
+    _meshNotDeleted = mesh.meshNotDeleted();
+    
     if(_isStructured)
     {
         _nxyz = mesh.getCellGridStructure() ;
@@ -1946,16 +1966,57 @@ Mesh::getBoundaryMesh ( void )  const
        return Mesh(_boundaryMesh);
 }
 
+Mesh 
+Mesh::getBoundaryGroupMesh ( std::string groupName, int nth_occurence )  const 
+{
+       //count occurences of groupName in known group name list
+       int count_occurences = std::count (_faceGroupNames.begin(),_faceGroupNames.end(),groupName);
+       
+    if( count_occurences ==0 )//No group found
+    {
+        cout<<"Mesh::getBoundaryGroupMesh Error : face group " << groupName << " does not exist"<<endl;
+        cout<<"Known face group names are " ;
+        for(int i=0; i<_faceGroupNames.size(); i++)
+                       cout<< _faceGroupNames[i]<<" ";
+               cout<< endl;
+        throw CdmathException("Required face group does not exist");
+    }
+    else if ( count_occurences <= nth_occurence)//Too many groups found
+    {
+        cout<<"Mesh::getBoundaryGroupMesh Error : "<<count_occurences<<" groups have name " << groupName<<", but you asked fo occurencer number "<<nth_occurence<<"which is too large"<<endl;
+        cout<<"Call function getBoundaryGroupMesh ( string groupName, int nth_group_match ) with nth_group_match between 0 and "<<count_occurences-1<<" to discriminate them "<<endl ;
+        throw CdmathException("Several face groups have the same name but you asked for an occurence that does not exsit");
+    }
+    else if( count_occurences >1 )//Wrning several groups found
+               cout<<"Warning : "<<count_occurences<<" groups have name " << groupName<<". Searching occurence number "<<nth_occurence<<endl;
+               
+       //search occurence of group name in known group name list
+       std::vector<std::string>::const_iterator it = _faceGroupNames.begin();
+       for (int i = 0; i<nth_occurence+1; i++)
+               it = std::find(it,_faceGroupNames.end(),groupName);
+               
+       return Mesh(_faceGroups[it - _faceGroupNames.begin()]); 
+}
+
 int 
 Mesh::getMaxNbNeighbours(EntityType type) const
 {
-    double result=0;
+    int result=0;
     
     if (type==CELLS)
        {
+               int nbNeib;//local number of neighbours
         for(int i=0; i<_numberOfCells; i++)
-            if(result < _cells[i].getNumberOfFaces())
-                result=_cells[i].getNumberOfFaces();
+        {
+            Cell Ci = _cells[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
+            
+            if(result < nbNeib)
+                result=nbNeib;
+               }
        }
     else if(type==NODES)
        {
@@ -1973,6 +2034,9 @@ void
 Mesh::writeVTK ( const std::string fileName ) const
 //----------------------------------------------------------------------
 {
+       if( !_isStructured && !_meshNotDeleted )
+               throw CdmathException("Mesh::writeVTK : Cannot save mesh : no MEDCouplingUMesh loaded");
+               
        string fname=fileName+".vtu";
        _mesh->writeVTK(fname.c_str()) ;
 }
@@ -1982,10 +2046,14 @@ void
 Mesh::writeMED ( const std::string fileName ) const
 //----------------------------------------------------------------------
 {
+       if( !_meshNotDeleted )
+               throw CdmathException("Mesh::writeMED : Cannot save mesh : no MEDCouplingUMesh loaded");
+               
        string fname=fileName+".med";
        MEDCouplingUMesh* mu=_mesh->buildUnstructured();
        MEDCoupling::WriteUMesh(fname.c_str(),mu,true);
 
+       /* Try to save mesh groups */
        //MEDFileUMesh meshMEDFile;
        //meshMEDFile.setMeshAtLevel(0,mu);
        //for(int i=0; i< _faceGroups.size(); i++)
@@ -1994,51 +2062,34 @@ Mesh::writeMED ( const std::string fileName ) const
        //MEDCoupling::meshMEDFile.write(fname.c_str(),2)       ;
        //else
        //MEDCoupling::meshMEDFile.write(fname.c_str(),1)       ;
-
-
-       mu->decrRef();
 }
 
 std::vector< int > 
-Mesh::getGroupFaceIds(std::string groupName) const
+Mesh::getFaceGroupIds(std::string groupName, bool isBoundaryGroup) const
 {
-    if( std::find(_faceGroupNames.begin(),_faceGroupNames.end(),groupName) == _faceGroupNames.end() )
+       std::vector<std::string>::const_iterator  it = std::find(_faceGroupNames.begin(),_faceGroupNames.end(),groupName);
+    if( it == _faceGroupNames.end() )
     {
         cout<<"Mesh::getGroupFaceIds Error : face group " << groupName << " does not exist"<<endl;
         throw CdmathException("Required face group does not exist");
     }
     else
-    {
-        std::vector< int > result(0);
-        for(int i=0; i<_boundaryFaceIds.size(); i++)
-        {
-            vector< std::string > groupNames = getFace(_boundaryFaceIds[i]).getGroupNames();
-            if( std::find(groupNames.begin(), groupNames.end(),groupName) != groupNames.end() )
-                result.push_back(_boundaryFaceIds[i]);
-        }
-        return result;
-    }
+    {   
+               return _faceGroupsIds[it-_faceGroupNames.begin()];
+       }
 }
 
 std::vector< int > 
-Mesh::getGroupNodeIds(std::string groupName) const
+Mesh::getNodeGroupIds(std::string groupName, bool isBoundaryGroup) const
 {
-    if( std::find(_nodeGroupNames.begin(),_nodeGroupNames.end(),groupName) == _nodeGroupNames.end() )
+       std::vector<std::string>::const_iterator  it = std::find(_nodeGroupNames.begin(),_nodeGroupNames.end(),groupName);
+    if( it == _nodeGroupNames.end() )
     {
         cout<<"Mesh::getGroupNodeIds Error : node group " << groupName << " does not exist"<<endl;
         throw CdmathException("Required node group does not exist");
     }
     else
-    {
-        std::vector< int > result(0);
-        for(int i=0; i<_boundaryFaceIds.size(); i++)
-        {
-            vector< std::string > groupNames = getNode(_boundaryFaceIds[i]).getGroupNames();
-            if( std::find(groupNames.begin(), groupNames.end(),groupName) != groupNames.end() )
-                result.push_back(_boundaryFaceIds[i]);
-        }
-        return result;
-    }
+        return _nodeGroupsIds[it-_nodeGroupNames.begin()];
 }
 
 void 
@@ -2046,6 +2097,10 @@ Mesh::setFaceGroupByIds(std::vector< int > faceIds, std::string groupName)
 {
     for(int i=0; i< faceIds.size(); i++)
         getFace(faceIds[i]).setGroupName(groupName);
+        
+    _faceGroupsIds.insert(  _faceGroupsIds.end(),faceIds);
+    _faceGroupNames.insert(_faceGroupNames.end(), groupName);
+    _faceGroups.insert(    _faceGroups.end(), NULL);
 }
 
 void 
@@ -2055,3 +2110,13 @@ Mesh::setNodeGroupByIds(std::vector< int > nodeIds, std::string groupName)
         getNode(nodeIds[i]).setGroupName(groupName);
 }
 
+void Mesh::deleteMEDCouplingUMesh()
+{ 
+       if(_meshNotDeleted) 
+       {
+               (_mesh.retn())->decrRef(); 
+               _meshNotDeleted=false;
+       } 
+       else 
+               throw CdmathException("Mesh::deleteMEDCouplingMesh() : mesh is not loaded");
+};