Salome HOME
Bug correction in the medcoupling <-> numpy layer
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingUMesh.cxx
index 697b876ab41c2ccf18c9ebb0565217f4e7f06206..0e1258aa705e26fb07fc97b83a8b73d6ff0ceaf0 100644 (file)
@@ -18,7 +18,7 @@
 //
 // Author : Anthony Geay (EDF R&D)
 
-#include "MEDCouplingUMesh.hxx"
+#include "MEDCouplingUMesh.txx"
 #include "MEDCouplingCMesh.hxx"
 #include "MEDCoupling1GTUMesh.hxx"
 #include "MEDCouplingFieldDouble.hxx"
@@ -77,7 +77,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::New(const std::string& meshName, int meshDim
  * Returns a new MEDCouplingUMesh which is a full copy of \a this one. No data is shared
  * between \a this and the new mesh.
  *  \return MEDCouplingUMesh * - a new instance of MEDCouplingMesh. The caller is to
- *          delete this mesh using decrRef() as it is no more needed. 
+ *          delete this mesh using decrRef() as it is no more needed.
  */
 MEDCouplingUMesh *MEDCouplingUMesh::deepCopy() const
 {
@@ -90,7 +90,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::deepCopy() const
  *  \param [in] recDeepCpy - if \a true, the copy is deep, else all data arrays of \a
  * this mesh are shared by the new mesh.
  *  \return MEDCouplingUMesh * - a new instance of MEDCouplingMesh. The caller is to
- *          delete this mesh using decrRef() as it is no more needed. 
+ *          delete this mesh using decrRef() as it is no more needed.
  */
 MEDCouplingUMesh *MEDCouplingUMesh::clone(bool recDeepCpy) const
 {
@@ -100,7 +100,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::clone(bool recDeepCpy) const
 /*!
  * This method behaves mostly like MEDCouplingUMesh::deepCopy method, except that only nodal connectivity arrays are deeply copied.
  * The coordinates are shared between \a this and the returned instance.
- * 
+ *
  * \return MEDCouplingUMesh * - A new object instance holding the copy of \a this (deep for connectivity, shallow for coordiantes)
  * \sa MEDCouplingUMesh::deepCopy
  */
@@ -301,9 +301,10 @@ void MEDCouplingUMesh::setMeshDimension(int meshDim)
 }
 
 /*!
- * Allocates memory to store an estimation of the given number of cells. The closer is the estimation to the number of cells effectively inserted,
- * the less will the library need to reallocate memory. If the number of cells to be inserted is not known simply put 0 to this parameter.
- * If a nodal connectivity previouly existed before the call of this method, it will be reset.
+ * Allocates memory to store an estimation of the given number of cells.
+ * The closer the estimation to the number of cells effectively inserted, the less need the library requires
+ * to reallocate memory. If the number of cells to be inserted is not known simply assign 0 to this parameter.
+ * If a nodal connectivity previously existed before the call of this method, it will be reset.
  *
  *  \param [in] nbOfCells - estimation of the number of cell \a this mesh will contain.
  *
@@ -339,7 +340,7 @@ void MEDCouplingUMesh::allocateCells(int nbOfCells)
  *  \param [in] type - type of cell to add.
  *  \param [in] size - number of nodes constituting this cell.
  *  \param [in] nodalConnOfCell - the connectivity of the cell to add.
- * 
+ *
  *  \if ENABLE_EXAMPLES
  *  \ref medcouplingcppexamplesUmeshStdBuild1 "Here is a C++ example".<br>
  *  \ref medcouplingpyexamplesUmeshStdBuild1 "Here is a Python example".
@@ -377,7 +378,7 @@ void MEDCouplingUMesh::insertNextCell(INTERP_KERNEL::NormalizedCellType type, in
 /*!
  * Compacts data arrays to release unused memory. This method is to be called after
  * finishing cell insertion using \a this->insertNextCell().
- * 
+ *
  *  \if ENABLE_EXAMPLES
  *  \ref medcouplingcppexamplesUmeshStdBuild1 "Here is a C++ example".<br>
  *  \ref medcouplingpyexamplesUmeshStdBuild1 "Here is a Python example".
@@ -403,7 +404,7 @@ MEDCouplingUMeshCellIterator *MEDCouplingUMesh::cellIterator()
 
 /*!
  * Entry point for iteration over cells groups geo types per geotypes. Warning the returned cell iterator should be deallocated.
- * If \a this is not so that that cells are grouped by geo types this method will throw an exception.
+ * If \a this is not so that the cells are grouped by geo types, this method will throw an exception.
  * In this case MEDCouplingUMesh::sortCellsInMEDFileFrmt or MEDCouplingUMesh::rearrange2ConsecutiveCellTypes methods for example can be called before invoking this method.
  * Useful for python users.
  */
@@ -431,7 +432,7 @@ std::set<INTERP_KERNEL::NormalizedCellType> MEDCouplingUMesh::getAllGeoTypes() c
  * having the same geometric type. So a same geometric type can appear more than once if the cells are not sorted per geometric type.
  *
  * \throw if connectivity in \a this is not correctly defined.
- *  
+ *
  * \sa MEDCouplingMesh::getAllGeoTypes
  */
 std::vector<INTERP_KERNEL::NormalizedCellType> MEDCouplingUMesh::getAllGeoTypesSorted() const
@@ -561,7 +562,7 @@ void MEDCouplingUMesh::checkFastEquivalWith(const MEDCouplingMesh *other, double
   MEDCouplingPointSet::checkFastEquivalWith(other,prec);
   const MEDCouplingUMesh *otherC=dynamic_cast<const MEDCouplingUMesh *>(other);
   if(!otherC)
-    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::checkFastEquivalWith : Two meshes are not not unstructured !"); 
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::checkFastEquivalWith : Two meshes are not not unstructured !");
 }
 
 /*!
@@ -573,15 +574,15 @@ void MEDCouplingUMesh::checkFastEquivalWith(const MEDCouplingMesh *other, double
  * \param [in,out] revNodalIndx - an array, of length \a this->getNumberOfNodes() + 1,
  *        dividing cell ids in \a revNodal into groups each referring to one
  *        node. Its every element (except the last one) is an index pointing to the
- *         first id of a group of cells. For example cells sharing the node #1 are 
- *        described by following range of indices: 
+ *         first id of a group of cells. For example cells sharing the node #1 are
+ *        described by following range of indices:
  *        [ \a revNodalIndx[1], \a revNodalIndx[2] ) and the cell ids are
  *        \a revNodal[ \a revNodalIndx[1] ], \a revNodal[ \a revNodalIndx[1] + 1], ...
  *        Number of cells sharing the *i*-th node is
  *        \a revNodalIndx[ *i*+1 ] - \a revNodalIndx[ *i* ].
  * \throw If the coordinates array is not set.
  * \throw If the nodal connectivity of cells is not defined.
- * 
+ *
  * \if ENABLE_EXAMPLES
  * \ref cpp_mcumesh_getReverseNodalConnectivity "Here is a C++ example".<br>
  * \ref  py_mcumesh_getReverseNodalConnectivity "Here is a Python example".
@@ -592,7 +593,7 @@ void MEDCouplingUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataA
   checkFullyDefined();
   int nbOfNodes(getNumberOfNodes());
   int *revNodalIndxPtr=(int *)malloc((nbOfNodes+1)*sizeof(int));
-  revNodalIndx->useArray(revNodalIndxPtr,true,C_DEALLOC,nbOfNodes+1,1);
+  revNodalIndx->useArray(revNodalIndxPtr,true,DeallocType::C_DEALLOC,nbOfNodes+1,1);
   std::fill(revNodalIndxPtr,revNodalIndxPtr+nbOfNodes+1,0);
   const int *conn(_nodal_connec->begin()),*connIndex(_nodal_connec_index->begin());
   int nbOfCells(getNumberOfCells()),nbOfEltsInRevNodal(0);
@@ -608,7 +609,7 @@ void MEDCouplingUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataA
     }
   std::transform(revNodalIndxPtr+1,revNodalIndxPtr+nbOfNodes+1,revNodalIndxPtr,revNodalIndxPtr+1,std::plus<int>());
   int *revNodalPtr=(int *)malloc((nbOfEltsInRevNodal)*sizeof(int));
-  revNodal->useArray(revNodalPtr,true,C_DEALLOC,nbOfEltsInRevNodal,1);
+  revNodal->useArray(revNodalPtr,true,DeallocType::C_DEALLOC,nbOfEltsInRevNodal,1);
   std::fill(revNodalPtr,revNodalPtr+nbOfEltsInRevNodal,-1);
   for(int eltId=0;eltId<nbOfCells;eltId++)
     {
@@ -627,7 +628,7 @@ void MEDCouplingUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataA
  * returned. The arrays \a desc and \a descIndx (\ref numbering-indirect) describe the descending connectivity,
  * i.e. enumerate cells of the result mesh bounding each cell of \a this mesh. The
  * arrays \a revDesc and \a revDescIndx (\ref numbering-indirect) describe the reverse descending connectivity,
- * i.e. enumerate cells of  \a this mesh bounded by each cell of the result mesh. 
+ * i.e. enumerate cells of  \a this mesh bounded by each cell of the result mesh.
  * \warning For speed reasons, this method does not check if node ids in the nodal
  *          connectivity correspond to the size of node coordinates array.
  * \warning Cells of the result mesh are \b not sorted by geometric type, hence,
@@ -657,7 +658,7 @@ void MEDCouplingUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataA
  *  \throw If the nodal connectivity of cells is node defined.
  *  \throw If \a desc == NULL || \a descIndx == NULL || \a revDesc == NULL || \a
  *         revDescIndx == NULL.
- * 
+ *
  *  \if ENABLE_EXAMPLES
  *  \ref cpp_mcumesh_buildDescendingConnectivity "Here is a C++ example".<br>
  *  \ref  py_mcumesh_buildDescendingConnectivity "Here is a Python example".
@@ -686,7 +687,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::explode3DMeshTo1D(DataArrayInt *desc, DataAr
 
 /*!
  * This method computes the micro edges constituting each cell in \a this. Micro edge is an edge for non quadratic cells. Micro edge is an half edge for quadratic cells.
- * This method works for both meshes with mesh dimenstion equal to 2 or 3. Dynamical cells are not supported (polygons, polyhedrons...)
+ * This method works for both meshes with mesh dimension equal to 2 or 3. Dynamical cells are not supported (polygons, polyhedrons...)
  * 
  * \sa explode3DMeshTo1D, buildDescendingConnectiviy
  */
@@ -717,7 +718,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::explodeMeshIntoMicroEdges(DataArrayInt *desc
  * cell with id #0 can't be negative, the array \a desc returns ids in FORTRAN mode,
  * i.e. cell ids are one-based.
  * Arrays \a revDesc and \a revDescIndx (\ref numbering-indirect) describe the reverse descending connectivity,
- * i.e. enumerate cells of  \a this mesh bounded by each cell of the result mesh. 
+ * i.e. enumerate cells of  \a this mesh bounded by each cell of the result mesh.
  * \warning For speed reasons, this method does not check if node ids in the nodal
  *          connectivity correspond to the size of node coordinates array.
  * \warning Cells of the result mesh are \b not sorted by geometric type, hence,
@@ -748,7 +749,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::explodeMeshIntoMicroEdges(DataArrayInt *desc
  *  \throw If the nodal connectivity of cells is node defined.
  *  \throw If \a desc == NULL || \a descIndx == NULL || \a revDesc == NULL || \a
  *         revDescIndx == NULL.
- * 
+ *
  *  \if ENABLE_EXAMPLES
  *  \ref cpp_mcumesh_buildDescendingConnectivity2 "Here is a C++ example".<br>
  *  \ref  py_mcumesh_buildDescendingConnectivity2 "Here is a Python example".
@@ -904,7 +905,7 @@ MCAuto<MEDCouplingUMesh> MEDCouplingUMesh::explodeIntoEdges(MCAuto<DataArrayInt>
  * The number of tuples is equal to the last values in \b neighborsIndx.
  * \param [out] neighborsIdx is an array of size this->getNumberOfCells()+1 newly allocated and should
  * be dealt by the caller. This arrays allow to use the first output parameter \b neighbors.
- * 
+ *
  * \sa MEDCouplingUMesh::computeEnlargedNeighborsOfNodes
  */
 void MEDCouplingUMesh::computeNeighborsOfNodes(DataArrayInt *&neighbors, DataArrayInt *&neighborsIdx) const
@@ -954,7 +955,7 @@ void MEDCouplingUMesh::computeNeighborsOfNodes(DataArrayInt *&neighbors, DataArr
  * Computes enlarged neighbors for each nodes in \a this. The behavior of this method is close to MEDCouplingUMesh::computeNeighborsOfNodes except that the neighborhood of each node is wider here.
  * A node j is considered to be in the neighborhood of i if and only if there is a cell in \a this containing in its nodal connectivity both i and j.
  * This method is useful to find ghost cells of a part of a mesh with a code based on fields on nodes.
- * 
+ *
  * \sa MEDCouplingUMesh::computeNeighborsOfNodes
  */
 void MEDCouplingUMesh::computeEnlargedNeighborsOfNodes(MCAuto<DataArrayInt> &neighbors, MCAuto<DataArrayInt>& neighborsIdx) const
@@ -975,7 +976,12 @@ void MEDCouplingUMesh::computeEnlargedNeighborsOfNodes(MCAuto<DataArrayInt> &nei
   {
     int *neighIdx(neighborsIdx->getPointer());
     for(std::vector< std::set<int> >::const_iterator it=st0.begin();it!=st0.end();it++,neighIdx++)
-      neighIdx[1]=neighIdx[0]+(*it).size()-1;
+      {
+        if ((*it).empty())
+          neighIdx[1]=neighIdx[0];
+        else
+          neighIdx[1]=neighIdx[0]+(*it).size()-1;
+      }
   }
   neighbors=DataArrayInt::New(); neighbors->alloc(neighborsIdx->back(),1);
   {
@@ -1288,9 +1294,9 @@ bool MEDCouplingUMesh::unPolyze()
 /*!
  * This method expects that spaceDimension is equal to 3 and meshDimension equal to 3.
  * This method performs operation only on polyhedrons in \b this. If no polyhedrons exists in \b this, \b this remains unchanged.
- * This method allows to merge if any coplanar 3DSurf cells that may appear in some polyhedrons cells. 
+ * This method allows to merge if any coplanar 3DSurf cells that may appear in some polyhedrons cells.
  *
- * \param [in] eps is a relative precision that allows to establish if some 3D plane are coplanar or not. This epsilon is used to recenter around origin to have maximal 
+ * \param [in] eps is a relative precision that allows to establish if some 3D plane are coplanar or not. This epsilon is used to recenter around origin to have maximal
  *             precision.
  */
 void MEDCouplingUMesh::simplifyPolyhedra(double eps)
@@ -1308,12 +1314,14 @@ void MEDCouplingUMesh::simplifyPolyhedra(double eps)
   connINew->alloc(nbOfCells+1,1);
   int *connINewPtr=connINew->getPointer(); *connINewPtr++=0;
   MCAuto<DataArrayInt> connNew=DataArrayInt::New(); connNew->alloc(0,1);
+  MCAuto<DataArrayInt> E_Fi(DataArrayInt::New()), E_F(DataArrayInt::New()), F_Ei(DataArrayInt::New()), F_E(DataArrayInt::New());
+  MCAuto<MEDCouplingUMesh> m_faces(buildDescendingConnectivity(E_F, E_Fi, F_E, F_Ei));
   bool changed=false;
   for(int i=0;i<nbOfCells;i++,connINewPtr++)
     {
       if(conn[index[i]]==(int)INTERP_KERNEL::NORM_POLYHED)
         {
-          SimplifyPolyhedronCell(eps,coords,conn+index[i],conn+index[i+1],connNew);
+          SimplifyPolyhedronCell(eps,coords, i,connNew, m_faces, E_Fi, E_F, F_Ei, F_E);
           changed=true;
         }
       else
@@ -1328,7 +1336,7 @@ void MEDCouplingUMesh::simplifyPolyhedra(double eps)
  * This method returns all node ids used in the connectivity of \b this. The data array returned has to be dealt by the caller.
  * The returned node ids are sorted ascendingly. This method is close to MEDCouplingUMesh::getNodeIdsInUse except
  * the format of the returned DataArrayInt instance.
- * 
+ *
  * \return a newly allocated DataArrayInt sorted ascendingly of fetched node ids.
  * \sa MEDCouplingUMesh::getNodeIdsInUse, areAllNodesFetched
  */
@@ -1378,12 +1386,12 @@ struct MEDCouplingAccVisit
 /*!
  * Finds nodes not used in any cell and returns an array giving a new id to every node
  * by excluding the unused nodes, for which the array holds -1. The result array is
- * a mapping in "Old to New" mode. 
+ * a mapping in "Old to New" mode.
  *  \param [out] nbrOfNodesInUse - number of node ids present in the nodal connectivity.
  *  \return DataArrayInt * - a new instance of DataArrayInt. Its length is \a
  *          this->getNumberOfNodes(). It holds for each node of \a this mesh either -1
  *          if the node is unused or a new id else. The caller is to delete this
- *          array using decrRef() as it is no more needed.  
+ *          array using decrRef() as it is no more needed.
  *  \throw If the coordinates array is not set.
  *  \throw If the nodal connectivity of cells is not defined.
  *  \throw If the nodal connectivity includes an invalid id.
@@ -1427,7 +1435,7 @@ DataArrayInt *MEDCouplingUMesh::getNodeIdsInUse(int& nbrOfNodesInUse) const
  * For each cell in \b this the number of nodes constituting cell is computed.
  * For each polyhedron cell, the sum of the number of nodes of each face constituting polyhedron cell is returned.
  * So for pohyhedrons some nodes can be counted several times in the returned result.
- * 
+ *
  * \return a newly allocated array
  * \sa MEDCouplingUMesh::computeEffectiveNbOfNodesPerCell
  */
@@ -1483,7 +1491,7 @@ DataArrayInt *MEDCouplingUMesh::computeEffectiveNbOfNodesPerCell() const
 /*!
  * This method returns a newly allocated array containing this->getNumberOfCells() tuples and 1 component.
  * For each cell in \b this the number of faces constituting (entity of dimension this->getMeshDimension()-1) cell is computed.
- * 
+ *
  * \return a newly allocated array
  */
 DataArrayInt *MEDCouplingUMesh::computeNbOfFacesPerCell() const
@@ -1506,10 +1514,10 @@ DataArrayInt *MEDCouplingUMesh::computeNbOfFacesPerCell() const
 /*!
  * Removes unused nodes (the node coordinates array is shorten) and returns an array
  * mapping between new and old node ids in "Old to New" mode. -1 values in the returned
- * array mean that the corresponding old node is no more used. 
+ * array mean that the corresponding old node is no more used.
  *  \return DataArrayInt * - a new instance of DataArrayInt of length \a
  *           this->getNumberOfNodes() before call of this method. The caller is to
- *           delete this array using decrRef() as it is no more needed. 
+ *           delete this array using decrRef() as it is no more needed.
  *  \throw If the coordinates array is not set.
  *  \throw If the nodal connectivity of cells is not defined.
  *  \throw If the nodal connectivity includes an invalid id.
@@ -1698,8 +1706,8 @@ int MEDCouplingUMesh::AreCellsEqualPolicy7(const int *conn, const int *connI, in
  * \param [in] startCellId specifies the cellId starting from which the equality computation will be carried out. By default it is 0, which it means that all cells in \a this will be scanned.
  * \param [out] commonCellsArr common cells ids (\ref numbering-indirect)
  * \param [out] commonCellsIArr common cells ids (\ref numbering-indirect)
- * \return the correspondance array old to new in a newly allocated array.
- * 
+ * \return the correspondence array old to new in a newly allocated array.
+ *
  */
 void MEDCouplingUMesh::findCommonCells(int compType, int startCellId, DataArrayInt *& commonCellsArr, DataArrayInt *& commonCellsIArr) const
 {
@@ -1795,7 +1803,7 @@ void MEDCouplingUMesh::FindCommonCellsAlg(int compType, int startCellId, const D
  * than \a this->getNumberOfCells() in the returned array means that there is no
  * corresponding cell in \a this mesh.
  * It is expected that \a this and \a other meshes share the same node coordinates
- * array, if it is not so an exception is thrown. 
+ * array, if it is not so an exception is thrown.
  *  \param [in] other - the mesh to compare with.
  *  \param [in] compType - specifies a cell comparison technique. For meaning of its
  *         valid values [0,1,2], see zipConnectivityTraducer().
@@ -1895,14 +1903,14 @@ MEDCouplingUMesh *MEDCouplingUMesh::mergeMyselfWithOnSameCoords(const MEDCouplin
 /*!
  * Build a sub part of \b this lying or not on the same coordinates than \b this (regarding value of \b keepCoords).
  * By default coordinates are kept. This method is close to MEDCouplingUMesh::buildPartOfMySelf except that here input
- * cellIds is not given explicitely but by a range python like.
- * 
+ * cellIds is not given explicitly but by a range python like.
+ *
  * \param start
  * \param end
  * \param step
  * \param keepCoords that specifies if you want or not to keep coords as this or zip it (see MEDCoupling::MEDCouplingUMesh::zipCoords). If true zipCoords is \b NOT called, if false, zipCoords is called.
  * \return a newly allocated
- * 
+ *
  * \warning This method modifies can generate an unstructured mesh whose cells are not sorted by geometric type order.
  * In view of the MED file writing, a renumbering of cells of returned unstructured mesh (using MEDCouplingUMesh::sortCellsInMEDFileFrmt) should be necessary.
  */
@@ -1935,7 +1943,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildPartOfMySelfSlice(int start, int end, i
  *         array of \a this mesh, else "free" nodes are removed from the result mesh
  *         by calling zipCoords().
  *  \return MEDCouplingUMesh * - a new instance of MEDCouplingUMesh. The caller is
- *         to delete this mesh using decrRef() as it is no more needed. 
+ *         to delete this mesh using decrRef() as it is no more needed.
  *  \throw If the coordinates array is not set.
  *  \throw If the nodal connectivity of cells is not defined.
  *  \throw If any cell id in the array \a begin is not valid.
@@ -2069,7 +2077,7 @@ void MEDCouplingUMesh::setPartOfMySelfSlice(int start, int end, int step, const
       MCAuto<DataArrayInt> arrOutAuto(arrOut),arrIOutAuto(arrIOut);
       setConnectivity(arrOut,arrIOut,true);
     }
-}                      
+}
 
 
 /*!
@@ -2079,14 +2087,14 @@ void MEDCouplingUMesh::setPartOfMySelfSlice(int start, int end, int step, const
  * specified node ids and the value of \a fullyIn parameter. If \a fullyIn ==\c true, a
  * cell is copied if its all nodes are in the array \a begin of node ids. If \a fullyIn
  * ==\c false, a cell is copied if any its node is in the array of node ids. The
- * created mesh shares the node coordinates array with \a this mesh. 
+ * created mesh shares the node coordinates array with \a this mesh.
  *  \param [in] begin - the array of node ids.
  *  \param [in] end - a pointer to the (last+1)-th element of \a begin.
  *  \param [in] fullyIn - if \c true, then cells whose all nodes are in the
  *         array \a begin are added, else cells whose any node is in the
  *         array \a begin are added.
  *  \return MEDCouplingUMesh * - new instance of MEDCouplingUMesh. The caller is
- *         to delete this mesh using decrRef() as it is no more needed. 
+ *         to delete this mesh using decrRef() as it is no more needed.
  *  \throw If the coordinates array is not set.
  *  \throw If the nodal connectivity of cells is not defined.
  *  \throw If any node id in \a begin is not valid.
@@ -2112,7 +2120,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildFacePartOfMySelfNode(const int *begin,
  *         array of \a this mesh, else "free" nodes are removed from the result mesh
  *         by calling zipCoords().
  *  \return MEDCouplingUMesh * - a new instance of MEDCouplingUMesh. The caller is
- *         to delete this mesh using decrRef() as it is no more needed. 
+ *         to delete this mesh using decrRef() as it is no more needed.
  *  \throw If the coordinates array is not set.
  *  \throw If the nodal connectivity of cells is not defined.
  *
@@ -2146,7 +2154,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildBoundaryMesh(bool keepCoords) const
 /*!
  * This method returns a newly created DataArrayInt instance containing ids of cells located in boundary.
  * A cell is detected to be on boundary if it contains one or more than one face having only one father.
- * This method makes the assumption that \a this is fully defined (coords,connectivity). If not an exception will be thrown. 
+ * This method makes the assumption that \a this is fully defined (coords,connectivity). If not an exception will be thrown.
  */
 DataArrayInt *MEDCouplingUMesh::findCellIdsOnBoundary() const
 {
@@ -2231,7 +2239,7 @@ void MEDCouplingUMesh::findCellIdsLyingOn(const MEDCouplingUMesh& otherDimM1OnSa
 /*!
  * This method computes the skin of \b this. That is to say the consituting meshdim-1 mesh is built and only the boundary subpart is
  * returned. This subpart of meshdim-1 mesh is built using meshdim-1 cells in it shared only one cell in \b this.
- * 
+ *
  * \return a newly allocated mesh lying on the same coordinates than \b this. The caller has to deal with returned mesh.
  */
 MEDCouplingUMesh *MEDCouplingUMesh::computeSkin() const
@@ -2277,7 +2285,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::buildUnstructured() const
  * This method expects that \b this and \b otherDimM1OnSameCoords share the same coordinates array.
  * otherDimM1OnSameCoords->getMeshDimension() is expected to be equal to this->getMeshDimension()-1.
  * This method searches for nodes needed to be duplicated. These nodes are nodes fetched by \b otherDimM1OnSameCoords which are not part of the boundary of \b otherDimM1OnSameCoords.
- * If a node is in the boundary of \b this \b and in the boundary of \b otherDimM1OnSameCoords this node is considerd as needed to be duplicated.
+ * If a node is in the boundary of \b this \b and in the boundary of \b otherDimM1OnSameCoords this node is considered as needed to be duplicated.
  * When the set of node ids \b nodeIdsToDuplicate is computed, cell ids in \b this is searched so that their connectivity includes at least 1 node in \b nodeIdsToDuplicate.
  *
  * \param [in] otherDimM1OnSameCoords a mesh lying on the same coords than \b this and with a mesh dimension equal to those of \b this minus 1. WARNING this input
@@ -2420,14 +2428,14 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On
 
 /*!
  * This method operates a modification of the connectivity and coords in \b this.
- * Every time that a node id in [ \b nodeIdsToDuplicateBg, \b nodeIdsToDuplicateEnd ) will append in nodal connectivity of \b this 
+ * Every time that a node id in [ \b nodeIdsToDuplicateBg, \b nodeIdsToDuplicateEnd ) will append in nodal connectivity of \b this
  * its ids will be modified to id this->getNumberOfNodes()+std::distance(nodeIdsToDuplicateBg,std::find(nodeIdsToDuplicateBg,nodeIdsToDuplicateEnd,id)).
- * More explicitely the renumber array in nodes is not explicitely given in old2new to avoid to build a big array of renumbering whereas typically few node ids needs to be
+ * More explicitly the renumber array in nodes is not explicitly given in old2new to avoid to build a big array of renumbering whereas typically few node ids needs to be
  * renumbered. The node id nodeIdsToDuplicateBg[0] will have id this->getNumberOfNodes()+0, node id nodeIdsToDuplicateBg[1] will have id this->getNumberOfNodes()+1,
  * node id nodeIdsToDuplicateBg[2] will have id this->getNumberOfNodes()+2...
- * 
+ *
  * As a consequence nodal connectivity array length will remain unchanged by this method, and nodal connectivity index array will remain unchanged by this method.
- * 
+ *
  * \param [in] nodeIdsToDuplicateBg begin of node ids (included) to be duplicated in connectivity only
  * \param [in] nodeIdsToDuplicateEnd end of node ids (excluded) to be duplicated in connectivity only
  */
@@ -2472,30 +2480,17 @@ void MEDCouplingUMesh::renumberNodesWithOffsetInConn(int offset)
  */
 void MEDCouplingUMesh::renumberNodesInConn(const INTERP_KERNEL::HashMap<int,int>& newNodeNumbersO2N)
 {
-  checkConnectivityFullyDefined();
-  int *conn(getNodalConnectivity()->getPointer());
-  const int *connIndex(getNodalConnectivityIndex()->getConstPointer());
-  int nbOfCells(getNumberOfCells());
-  for(int i=0;i<nbOfCells;i++)
-    for(int iconn=connIndex[i]+1;iconn!=connIndex[i+1];iconn++)
-      {
-        int& node=conn[iconn];
-        if(node>=0)//avoid polyhedron separator
-          {
-            INTERP_KERNEL::HashMap<int,int>::const_iterator it(newNodeNumbersO2N.find(node));
-            if(it!=newNodeNumbersO2N.end())
-              {
-                node=(*it).second;
-              }
-            else
-              {
-                std::ostringstream oss; oss << "MEDCouplingUMesh::renumberNodesInConn(map) : presence in connectivity for cell #" << i << " of node #" << node << " : Not in map !";
-                throw INTERP_KERNEL::Exception(oss.str());
-              }
-          }
-      }
-  _nodal_connec->declareAsNew();
-  updateTime();
+  this->renumberNodesInConnT< INTERP_KERNEL::HashMap<int,int> >(newNodeNumbersO2N);
+}
+
+/*!
+ *  Same than renumberNodesInConn(const int *) except that here the format of old-to-new traducer is using map instead
+ *  of array. This method is dedicated for renumbering from a big set of nodes the a tiny set of nodes which is the case during extraction
+ *  of a big mesh.
+ */
+void MEDCouplingUMesh::renumberNodesInConn(const std::map<int,int>& newNodeNumbersO2N)
+{
+  this->renumberNodesInConnT< std::map<int,int> >(newNodeNumbersO2N);
 }
 
 /*!
@@ -2504,7 +2499,7 @@ void MEDCouplingUMesh::renumberNodesInConn(const INTERP_KERNEL::HashMap<int,int>
  * This method is a generalization of shiftNodeNumbersInConn().
  *  \warning This method performs no check of validity of new ids. **Use it with care !**
  *  \param [in] newNodeNumbersO2N - a permutation array, of length \a
- *         this->getNumberOfNodes(), in "Old to New" mode. 
+ *         this->getNumberOfNodes(), in "Old to New" mode.
  *         See \ref numbering for more info on renumbering modes.
  *  \throw If the nodal connectivity of cells is not defined.
  *
@@ -2536,7 +2531,7 @@ void MEDCouplingUMesh::renumberNodesInConn(const int *newNodeNumbersO2N)
  * This method renumbers nodes \b in \b connectivity \b only \b without \b any \b reference \b to \b coords.
  * This method performs no check on the fact that new coordinate ids are valid. \b Use \b it \b with \b care !
  * This method is an specialization of \ref MEDCoupling::MEDCouplingUMesh::renumberNodesInConn "renumberNodesInConn method".
- * 
+ *
  * \param [in] delta specifies the shift size applied to nodeId in nodal connectivity in \b this.
  */
 void MEDCouplingUMesh::shiftNodeNumbersInConn(int delta)
@@ -2561,18 +2556,18 @@ void MEDCouplingUMesh::shiftNodeNumbersInConn(int delta)
 /*!
  * This method operates a modification of the connectivity in \b this.
  * Coordinates are \b NOT considered here and will remain unchanged by this method. this->_coords can ever been null for the needs of this method.
- * Every time that a node id in [ \b nodeIdsToDuplicateBg, \b nodeIdsToDuplicateEnd ) will append in nodal connectivity of \b this 
+ * Every time that a node id in [ \b nodeIdsToDuplicateBg, \b nodeIdsToDuplicateEnd ) will append in nodal connectivity of \b this
  * its ids will be modified to id offset+std::distance(nodeIdsToDuplicateBg,std::find(nodeIdsToDuplicateBg,nodeIdsToDuplicateEnd,id)).
- * More explicitely the renumber array in nodes is not explicitely given in old2new to avoid to build a big array of renumbering whereas typically few node ids needs to be
+ * More explicitly the renumber array in nodes is not explicitly given in old2new to avoid to build a big array of renumbering whereas typically few node ids needs to be
  * renumbered. The node id nodeIdsToDuplicateBg[0] will have id offset+0, node id nodeIdsToDuplicateBg[1] will have id offset+1,
  * node id nodeIdsToDuplicateBg[2] will have id offset+2...
- * 
+ *
  * As a consequence nodal connectivity array length will remain unchanged by this method, and nodal connectivity index array will remain unchanged by this method.
  * As an another consequense after the call of this method \b this can be transiently non cohrent.
- * 
+ *
  * \param [in] nodeIdsToDuplicateBg begin of node ids (included) to be duplicated in connectivity only
  * \param [in] nodeIdsToDuplicateEnd end of node ids (excluded) to be duplicated in connectivity only
- * \param [in] offset the offset applied to all node ids in connectivity that are in [ \a nodeIdsToDuplicateBg, \a nodeIdsToDuplicateEnd ). 
+ * \param [in] offset the offset applied to all node ids in connectivity that are in [ \a nodeIdsToDuplicateBg, \a nodeIdsToDuplicateEnd ).
  */
 void MEDCouplingUMesh::duplicateNodesInConn(const int *nodeIdsToDuplicateBg, const int *nodeIdsToDuplicateEnd, int offset)
 {
@@ -2611,7 +2606,7 @@ void MEDCouplingUMesh::duplicateNodesInConn(const int *nodeIdsToDuplicateBg, con
  * If 'check' equals false the method will not check the content of [ \a old2NewBg ; \a old2NewEnd ).
  * To avoid any throw of SIGSEGV when 'check' equals false, the elements in [ \a old2NewBg ; \a old2NewEnd ) should be unique and
  * should be contained in[0;this->getNumberOfCells()).
- * 
+ *
  * \param [in] old2NewBg is expected to be a dynamically allocated pointer of size at least equal to this->getNumberOfCells()
  * \param check
  */
@@ -2625,7 +2620,7 @@ void MEDCouplingUMesh::renumberCells(const int *old2NewBg, bool check)
   //
   const int *conn=_nodal_connec->getConstPointer();
   const int *connI=_nodal_connec_index->getConstPointer();
-  MCAuto<DataArrayInt> o2n=DataArrayInt::New(); o2n->useArray(array,false,C_DEALLOC,nbCells,1);
+  MCAuto<DataArrayInt> o2n=DataArrayInt::New(); o2n->useArray(array,false,DeallocType::C_DEALLOC,nbCells,1);
   MCAuto<DataArrayInt> n2o=o2n->invertArrayO2N2N2O(nbCells);
   const int *n2oPtr=n2o->begin();
   MCAuto<DataArrayInt> newConn=DataArrayInt::New();
@@ -2657,13 +2652,13 @@ void MEDCouplingUMesh::renumberCells(const int *old2NewBg, bool check)
  * Finds cells whose bounding boxes intersect a given bounding box.
  *  \param [in] bbox - an array defining the bounding box via coordinates of its
  *         extremum points in "no interlace" mode, i.e. xMin, xMax, yMin, yMax, zMin,
- *         zMax (if in 3D). 
+ *         zMax (if in 3D).
  *  \param [in] eps - a factor used to increase size of the bounding box of cell
  *         before comparing it with \a bbox. This factor is multiplied by the maximal
  *         extent of the bounding box of cell to produce an addition to this bounding box.
  *  \return DataArrayInt * - a new instance of DataArrayInt holding ids for found
  *         cells. The caller is to delete this array using decrRef() as it is no more
- *         needed. 
+ *         needed.
  *  \throw If the coordinates array is not set.
  *  \throw If the nodal connectivity of cells is not defined.
  *
@@ -2934,11 +2929,11 @@ std::string MEDCouplingUMesh::reprConnectivityOfThis() const
 }
 
 /*!
- * This method builds a newly allocated instance (with the same name than \a this) that the caller has the responsability to deal with.
+ * This method builds a newly allocated instance (with the same name than \a this) that the caller has the responsibility to deal with.
  * This method returns an instance with all arrays allocated (connectivity, connectivity index, coordinates)
  * but with length of these arrays set to 0. It allows to define an "empty" mesh (with nor cells nor nodes but compliant with
  * some algos).
- * 
+ *
  * This method expects that \a this has a mesh dimension set and higher or equal to 0. If not an exception will be thrown.
  * This method analyzes the 3 arrays of \a this. For each the following behaviour is done : if the array is null a newly one is created
  * with number of tuples set to 0, if not the array is taken as this in the returned instance.
@@ -2996,11 +2991,11 @@ int MEDCouplingUMesh::getNumberOfNodesInCell(int cellId) const
 
 /*!
  * Returns types of cells of the specified part of \a this mesh.
- * This method avoids computing sub-mesh explicitely to get its types.
+ * This method avoids computing sub-mesh explicitly to get its types.
  *  \param [in] begin - an array of cell ids of interest.
  *  \param [in] end - the end of \a begin, i.e. a pointer to its (last+1)-th element.
  *  \return std::set<INTERP_KERNEL::NormalizedCellType> - a set of enumeration items
- *         describing the cell types. 
+ *         describing the cell types.
  *  \throw If the coordinates array is not set.
  *  \throw If the nodal connectivity of cells is not defined.
  *  \sa getAllGeoTypes()
@@ -3019,10 +3014,10 @@ std::set<INTERP_KERNEL::NormalizedCellType> MEDCouplingUMesh::getTypesOfPart(con
 /*!
  * Defines the nodal connectivity using given connectivity arrays in \ref numbering-indirect format.
  * Optionally updates
- * a set of types of cells constituting \a this mesh. 
+ * a set of types of cells constituting \a this mesh.
  * This method is for advanced users having prepared their connectivity before. For
  * more info on using this method see \ref MEDCouplingUMeshAdvBuild.
- *  \param [in] conn - the nodal connectivity array. 
+ *  \param [in] conn - the nodal connectivity array.
  *  \param [in] connIndex - the nodal connectivity index array.
  *  \param [in] isComputingTypes - if \c true, the set of types constituting \a this
  *         mesh is updated.
@@ -3069,12 +3064,12 @@ void MEDCouplingUMesh::computeTypes()
 
 
 /*!
- * Returns a number of cells constituting \a this mesh. 
+ * Returns a number of cells constituting \a this mesh.
  *  \return int - the number of cells in \a this mesh.
  *  \throw If the nodal connectivity of cells is not defined.
  */
 std::size_t MEDCouplingUMesh::getNumberOfCells() const
-{ 
+{
   if(_nodal_connec_index)
     return _nodal_connec_index->getNumberOfTuples()-1;
   else
@@ -3250,7 +3245,7 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureField(bool isAbs) const
  *  \param [in] end - the end of \a begin, i.e. a pointer to its (last+1)-th element.
  *  \return DataArrayDouble * - a new instance of DataArrayDouble. The caller is to
  *          delete this array using decrRef() as it is no more needed.
- * 
+ *
  *  \if ENABLE_EXAMPLES
  *  \ref cpp_mcumesh_getPartMeasureField "Here is a C++ example".<br>
  *  \ref  py_mcumesh_getPartMeasureField "Here is a Python example".
@@ -3336,10 +3331,10 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getMeasureFieldOnNode(bool isAbs) cons
  * mesh. The returned normal vectors to each cell have a norm2 equal to 1.
  * The computed vectors have <em> this->getMeshDimension()+1 </em> components
  * and are normalized.
- * <br> \a this can be either 
- * - a  2D mesh in 2D or 3D space or 
+ * <br> \a this can be either
+ * - a  2D mesh in 2D or 3D space or
  * - an 1D mesh in 2D space.
- * 
+ *
  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble on
  *          cells and one time. The caller is to delete this field using decrRef() as
  *          it is no more needed.
@@ -3406,10 +3401,10 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::buildOrthogonalField() const
  * Returns a new MEDCouplingFieldDouble holding normal vectors to specified cells of
  * \a this mesh. The computed vectors have <em> this->getMeshDimension()+1 </em> components
  * and are normalized.
- * <br> \a this can be either 
- * - a  2D mesh in 2D or 3D space or 
+ * <br> \a this can be either
+ * - a  2D mesh in 2D or 3D space or
  * - an 1D mesh in 2D space.
- * 
+ *
  * This method avoids building explicitly a part of \a this mesh to perform the work.
  *  \param [in] begin - an array of cell ids of interest.
  *  \param [in] end - the end of \a begin, i.e. a pointer to its (last+1)-th element.
@@ -3533,7 +3528,7 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::buildDirectionVectorField() const
  *         using decrRef() as it is no more needed.
  *  \return MEDCouplingUMesh * - a new instance of MEDCouplingUMesh. This mesh does
  *         not share the node coordinates array with \a this mesh. The caller is to
- *         delete this mesh using decrRef() as it is no more needed.  
+ *         delete this mesh using decrRef() as it is no more needed.
  *  \throw If the coordinates array is not set.
  *  \throw If the nodal connectivity of cells is not defined.
  *  \throw If \a this->getMeshDimension() != 3 or \a this->getSpaceDimension() != 3.
@@ -3599,7 +3594,7 @@ the result mesh.
  *  \return MEDCouplingUMesh * - a new instance of MEDCouplingUMesh. This is an 1D
  *         mesh in 3D space. This mesh does not share the node coordinates array with
  *         \a this mesh. The caller is to delete this mesh using decrRef() as it is
- *         no more needed. 
+ *         no more needed.
  *  \throw If the coordinates array is not set.
  *  \throw If the nodal connectivity of cells is not defined.
  *  \throw If \a this->getMeshDimension() != 2 or \a this->getSpaceDimension() != 3.
@@ -3839,7 +3834,7 @@ DataArrayInt *MEDCouplingUMesh::getCellIdsCrossingPlane(const double *origin, co
  * If not an exception will thrown. If this is an empty mesh with no cell an exception will be thrown too.
  * No consideration of coordinate is done by this method.
  * A 1D mesh is said contiguous if : a cell i with nodal connectivity (k,p) the cell i+1 the nodal connectivity should be (p,m)
- * If not false is returned. In case that false is returned a call to MEDCoupling::MEDCouplingUMesh::mergeNodes could be usefull.
+ * If not false is returned. In case that false is returned a call to MEDCoupling::MEDCouplingUMesh::mergeNodes could be useful.
  */
 bool MEDCouplingUMesh::isContiguous1D() const
 {
@@ -3899,7 +3894,7 @@ void MEDCouplingUMesh::project1D(const double *pt, const double *v, double eps,
 }
 
 /*!
- * This method computes the distance from a point \a pt to \a this and the first \a cellId in \a this corresponding to the returned distance. 
+ * This method computes the distance from a point \a pt to \a this and the first \a cellId in \a this corresponding to the returned distance.
  * \a this is expected to be a mesh so that its space dimension is equal to its
  * mesh dimension + 1. Furthermore only mesh dimension 1 and 2 are supported for the moment.
  * Distance from \a ptBg to \a ptEnd is expected to be equal to the space dimension. \a this is also expected to be fully defined (connectivity and coordinates).
@@ -3930,7 +3925,7 @@ double MEDCouplingUMesh::distanceToPoint(const double *ptBg, const double *ptEnd
   if((int)std::distance(ptBg,ptEnd)!=spaceDim)
     { std::ostringstream oss; oss << "MEDCouplingUMesh::distanceToPoint : input point has to have dimension equal to the space dimension of this (" << spaceDim << ") !"; throw INTERP_KERNEL::Exception(oss.str()); }
   DataArrayInt *ret1=0;
-  MCAuto<DataArrayDouble> pts=DataArrayDouble::New(); pts->useArray(ptBg,false,C_DEALLOC,1,spaceDim);
+  MCAuto<DataArrayDouble> pts=DataArrayDouble::New(); pts->useArray(ptBg,false,DeallocType::C_DEALLOC,1,spaceDim);
   MCAuto<DataArrayDouble> ret0=distanceToPoints(pts,ret1);
   MCAuto<DataArrayInt> ret1Safe(ret1);
   cellId=*ret1Safe->begin();
@@ -3939,11 +3934,11 @@ double MEDCouplingUMesh::distanceToPoint(const double *ptBg, const double *ptEnd
 
 /*!
  * This method computes the distance from each point of points serie \a pts (stored in a DataArrayDouble in which each tuple represents a point)
- *  to \a this  and the first \a cellId in \a this corresponding to the returned distance. 
+ *  to \a this  and the first \a cellId in \a this corresponding to the returned distance.
  * WARNING, if there is some orphan nodes in \a this (nodes not fetched by any cells in \a this ( see MEDCouplingUMesh::zipCoords ) ) these nodes will ** not ** been taken
  * into account in this method. Only cells and nodes lying on them are considered in the algorithm (even if one of these orphan nodes is closer than returned distance).
  * A user that needs to consider orphan nodes should invoke DataArrayDouble::minimalDistanceTo method on the coordinates array of \a this.
- * 
+ *
  * \a this is expected to be a mesh so that its space dimension is equal to its
  * mesh dimension + 1. Furthermore only mesh dimension 1 and 2 are supported for the moment.
  * Number of components of \a pts is expected to be equal to the space dimension. \a this is also expected to be fully defined (connectivity and coordinates).
@@ -4024,13 +4019,13 @@ DataArrayDouble *MEDCouplingUMesh::distanceToPoints(const DataArrayDouble *pts,
 /// @endcond
 
 /*!
- * Finds cells in contact with a ball (i.e. a point with precision). 
+ * Finds cells in contact with a ball (i.e. a point with precision).
  * For speed reasons, the INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6 and INTERP_KERNEL::NORM_QUAD8 cells are considered as convex cells to detect if a point is IN or OUT.
  * If it is not the case, please change their types to INTERP_KERNEL::NORM_POLYGON or INTERP_KERNEL::NORM_QPOLYG before invoking this method.
  *
  * \warning This method is suitable if the caller intends to evaluate only one
  *          point, for more points getCellsContainingPoints() is recommended as it is
- *          faster. 
+ *          faster.
  *  \param [in] pos - array of coordinates of the ball central point.
  *  \param [in] eps - ball radius.
  *  \return int - a smallest id of cells being in contact with the ball, -1 in case
@@ -4053,7 +4048,7 @@ int MEDCouplingUMesh::getCellContainingPoint(const double *pos, double eps) cons
  * If it is not the case, please change their types to INTERP_KERNEL::NORM_POLYGON or INTERP_KERNEL::NORM_QPOLYG before invoking this method.
  * \warning This method is suitable if the caller intends to evaluate only one
  *          point, for more points getCellsContainingPoints() is recommended as it is
- *          faster. 
+ *          faster.
  *  \param [in] pos - array of coordinates of the ball central point.
  *  \param [in] eps - ball radius.
  *  \param [out] elts - vector returning ids of the found cells. It is cleared
@@ -4073,6 +4068,45 @@ void MEDCouplingUMesh::getCellsContainingPoint(const double *pos, double eps, st
   elts.clear(); elts.insert(elts.end(),eltsUg->begin(),eltsUg->end());
 }
 
+void MEDCouplingUMesh::getCellsContainingPointsZeAlg(const double *pos, int nbOfPoints, double eps,
+                                                     MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& eltsIndex,
+                                                     std::function<bool(INTERP_KERNEL::NormalizedCellType,int)> sensibilityTo2DQuadraticLinearCellsFunc) const
+{
+  int spaceDim(getSpaceDimension()),mDim(getMeshDimension());
+  if(spaceDim==3)
+    {
+      if(mDim==3)
+        {
+          const double *coords=_coords->getConstPointer();
+          getCellsContainingPointsAlg<3>(coords,pos,nbOfPoints,eps,elts,eltsIndex,sensibilityTo2DQuadraticLinearCellsFunc);
+        }
+      else
+        throw INTERP_KERNEL::Exception("For spaceDim==3 only meshDim==3 implemented for getelementscontainingpoints !");
+    }
+  else if(spaceDim==2)
+    {
+      if(mDim==2)
+        {
+          const double *coords=_coords->getConstPointer();
+          getCellsContainingPointsAlg<2>(coords,pos,nbOfPoints,eps,elts,eltsIndex,sensibilityTo2DQuadraticLinearCellsFunc);
+        }
+      else
+        throw INTERP_KERNEL::Exception("For spaceDim==2 only meshDim==2 implemented for getelementscontainingpoints !");
+    }
+  else if(spaceDim==1)
+    {
+      if(mDim==1)
+        {
+          const double *coords=_coords->getConstPointer();
+          getCellsContainingPointsAlg<1>(coords,pos,nbOfPoints,eps,elts,eltsIndex,sensibilityTo2DQuadraticLinearCellsFunc);
+        }
+      else
+        throw INTERP_KERNEL::Exception("For spaceDim==1 only meshDim==1 implemented for getelementscontainingpoints !");
+    }
+  else
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getCellsContainingPoints : not managed for mdim not in [1,2,3] !");
+}
+
 /*!
  * Finds cells in contact with several balls (i.e. points with precision).
  * This method is an extension of getCellContainingPoint() and
@@ -4081,7 +4115,7 @@ void MEDCouplingUMesh::getCellsContainingPoint(const double *pos, double eps, st
  * If it is not the case, please change their types to INTERP_KERNEL::NORM_POLYGON or INTERP_KERNEL::NORM_QPOLYG before invoking this method.
  *  \param [in] pos - an array of coordinates of points in full interlace mode :
  *         X0,Y0,Z0,X1,Y1,Z1,... Size of the array must be \a
- *         this->getSpaceDimension() * \a nbOfPoints 
+ *         this->getSpaceDimension() * \a nbOfPoints
  *  \param [in] nbOfPoints - number of points to locate within \a this mesh.
  *  \param [in] eps - radius of balls (i.e. the precision).
  *  \param [out] elts - vector returning ids of found cells.
@@ -4105,44 +4139,21 @@ void MEDCouplingUMesh::getCellsContainingPoint(const double *pos, double eps, st
 void MEDCouplingUMesh::getCellsContainingPoints(const double *pos, int nbOfPoints, double eps,
                                                 MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& eltsIndex) const
 {
-  int spaceDim=getSpaceDimension();
-  int mDim=getMeshDimension();
-  if(spaceDim==3)
-    {
-      if(mDim==3)
-        {
-          const double *coords=_coords->getConstPointer();
-          getCellsContainingPointsAlg<3>(coords,pos,nbOfPoints,eps,elts,eltsIndex);
-        }
-      /*else if(mDim==2)
-        {
+  auto yesImSensibleTo2DQuadraticLinearCellsFunc([](INTERP_KERNEL::NormalizedCellType ct, int mdim) { return INTERP_KERNEL::CellModel::GetCellModel(ct).isQuadratic() && mdim == 2; } );
+  this->getCellsContainingPointsZeAlg(pos,nbOfPoints,eps,elts,eltsIndex,yesImSensibleTo2DQuadraticLinearCellsFunc);
+}
 
-        }*/
-      else
-        throw INTERP_KERNEL::Exception("For spaceDim==3 only meshDim==3 implemented for getelementscontainingpoints !");
-    }
-  else if(spaceDim==2)
-    {
-      if(mDim==2)
-        {
-          const double *coords=_coords->getConstPointer();
-          getCellsContainingPointsAlg<2>(coords,pos,nbOfPoints,eps,elts,eltsIndex);
-        }
-      else
-        throw INTERP_KERNEL::Exception("For spaceDim==2 only meshDim==2 implemented for getelementscontainingpoints !");
-    }
-  else if(spaceDim==1)
-    {
-      if(mDim==1)
-        {
-          const double *coords=_coords->getConstPointer();
-          getCellsContainingPointsAlg<1>(coords,pos,nbOfPoints,eps,elts,eltsIndex);
-        }
-      else
-        throw INTERP_KERNEL::Exception("For spaceDim==1 only meshDim==1 implemented for getelementscontainingpoints !");
-    }
-  else
-    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getCellsContainingPoints : not managed for mdim not in [1,2,3] !");
+/*!
+ * Behaves like MEDCouplingMesh::getCellsContainingPoints for cells in \a this that are linear.
+ * For quadratic cells in \a this, this method behaves by just considering linear part of cells.
+ * This method is here only for backward compatibility (interpolation GaussPoints to GaussPoints).
+ * 
+ * \sa MEDCouplingUMesh::getCellsContainingPoints, MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss
+ */
+void MEDCouplingUMesh::getCellsContainingPointsLinearPartOnlyOnNonDynType(const double *pos, int nbOfPoints, double eps, MCAuto<DataArrayInt>& elts, MCAuto<DataArrayInt>& eltsIndex) const
+{
+  auto noImNotSensibleTo2DQuadraticLinearCellsFunc([](INTERP_KERNEL::NormalizedCellType,int) { return false; } );
+  this->getCellsContainingPointsZeAlg(pos,nbOfPoints,eps,elts,eltsIndex,noImNotSensibleTo2DQuadraticLinearCellsFunc);
 }
 
 /*!
@@ -4186,7 +4197,7 @@ void MEDCouplingUMesh::checkButterflyCells(std::vector<int>& cells, double eps)
  *
  * This method expects that space dimension is equal to 2 and mesh dimension is equal to 2 too. If it is not the case an INTERP_KERNEL::Exception will be thrown.
  * This method works only for linear 2D cells. If there is any of non linear cells (INTERP_KERNEL::NORM_QUAD8 for example) an INTERP_KERNEL::Exception will be thrown too.
- * 
+ *
  * For each 2D linear cell in \b this, this method builds the convex envelop (or the convex hull) of the current cell.
  * This convex envelop is computed using Jarvis march algorithm.
  * The coordinates and the number of cells of \b this remain unchanged on invocation of this method.
@@ -4239,7 +4250,7 @@ DataArrayInt *MEDCouplingUMesh::convexEnvelop2D()
  *   the 3 preceding points of the 1D mesh. The center of the arc is the center of rotation for each level, the rotation is done
  *   along an axis normal to the plane containing the arc, and finally the angle of rotation is defined by the first two points on the
  *   arc.
- * \return an unstructured mesh with meshDim==3 and spaceDim==3. The returned mesh has the same coords than \a this.  
+ * \return an unstructured mesh with meshDim==3 and spaceDim==3. The returned mesh has the same coords than \a this.
  */
 MEDCouplingUMesh *MEDCouplingUMesh::buildExtrudedMesh(const MEDCouplingUMesh *mesh1D, int policy)
 {
@@ -4393,11 +4404,11 @@ void MEDCouplingUMesh::convertQuadraticCellsToLinear()
  * or to INTERP_KERNEL::NORM_TRI7 if \a conversionType is equal to 1. All non linear cells and polyhedron in \a this are let untouched.
  * Contrary to MEDCouplingUMesh::convertQuadraticCellsToLinear method, the coordinates in \a this can be become bigger. All created nodes will be put at the
  * end of the existing coordinates.
- * 
+ *
  * \param [in] conversionType specifies the type of conversion expected. Only 0 (default) and 1 are supported presently. 0 those that creates the 'most' simple
  *             corresponding quadratic cells. 1 is those creating the 'most' complex.
  * \return a newly created DataArrayInt instance that the caller should deal with containing cell ids of converted cells.
- * 
+ *
  * \throw if \a this is not fully defined. It throws too if \a conversionType is not in [0,1].
  *
  * \sa MEDCouplingUMesh::convertQuadraticCellsToLinear
@@ -4491,16 +4502,6 @@ void MEDCouplingUMesh::tessellate2D(double eps)
       throw INTERP_KERNEL::Exception("MEDCouplingUMesh::tessellate2D : mesh dimension must be in [1,2] !");
     }
 }
-/*!
- * Tessellates \a this 1D mesh in 2D space by dividing not straight quadratic edges.
- * \warning This method can lead to a huge amount of nodes if \a eps is very low.
- *  \param [in] eps - specifies the maximal angle (in radian) between 2 sub-edges of
- *         a sub-divided edge.
- *  \throw If the coordinates array is not set.
- *  \throw If the nodal connectivity of cells is not defined.
- *  \throw If \a this->getMeshDimension() != 1.
- *  \throw If \a this->getSpaceDimension() != 2.
- */
 
 #if 0
 /*!
@@ -4562,7 +4563,7 @@ void MEDCouplingUMesh::splitSomeEdgesOf2DMesh(const DataArrayInt *nodeIdsToAdd,
  *
  *  \throw If \a policy is 0 or 1 and \a this->getMeshDimension() != 2.
  *  \throw If \a policy is INTERP_KERNEL::PLANAR_FACE_5 or INTERP_KERNEL::PLANAR_FACE_6
- *          and \a this->getMeshDimension() != 3. 
+ *          and \a this->getMeshDimension() != 3.
  *  \throw If \a policy is not one of the four discussed above.
  *  \throw If the nodal connectivity of cells is not defined.
  * \sa MEDCouplingUMesh::tetrahedrize, MEDCoupling1SGTUMesh::sortHexa8EachOther
@@ -4621,8 +4622,9 @@ bool MEDCouplingUMesh::areOnlySimplexCells() const
 /*!
  * Converts degenerated 2D or 3D linear cells of \a this mesh into cells of simpler
  * type. For example an INTERP_KERNEL::NORM_QUAD4 cell having only three unique nodes in
- * its connectivity is transformed into an INTERP_KERNEL::NORM_TRI3 cell. This method
- * does \b not perform geometrical checks and checks only nodal connectivity of cells,
+ * its connectivity is transformed into an INTERP_KERNEL::NORM_TRI3 cell.
+ * Quadratic cells in 2D are also handled. In those cells edges where start=end=midpoint are removed.
+ * This method does \b not perform geometrical checks and checks only nodal connectivity of cells,
  * so it can be useful to call mergeNodes() before calling this method.
  *  \throw If \a this->getMeshDimension() <= 1.
  *  \throw If the coordinates array is not set.
@@ -4659,12 +4661,114 @@ void MEDCouplingUMesh::convertDegeneratedCells()
   computeTypes();
 }
 
+/*!
+ * Same as MEDCouplingUMesh::convertDegeneratedCells() plus deletion of the flat cells.
+ * A cell is flat in the following cases:
+ *   - for a linear cell, all points in the connectivity are equal
+ *   - for a quadratic cell, either the above, or a quadratic polygon with two (linear) points and two
+ *   identical quadratic points
+ * \return a new instance of DataArrayInt holding ids of removed cells. The caller is to delete
+ *      this array using decrRef() as it is no more needed.
+ */
+DataArrayInt *MEDCouplingUMesh::convertDegeneratedCellsAndRemoveFlatOnes()
+{
+  checkFullyDefined();
+  if(getMeshDimension()<=1)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::convertDegeneratedCells works on umeshes with meshdim equals to 2 or 3 !");
+  int nbOfCells=getNumberOfCells();
+  MCAuto<DataArrayInt> ret(DataArrayInt::New()); ret->alloc(0,1);
+  if(nbOfCells<1)
+    return ret.retn();
+  int initMeshLgth=getNodalConnectivityArrayLen();
+  int *conn=_nodal_connec->getPointer();
+  int *index=_nodal_connec_index->getPointer();
+  int posOfCurCell=0;
+  int newPos=0;
+  int lgthOfCurCell, nbDelCells(0);
+  for(int i=0;i<nbOfCells;i++)
+    {
+      lgthOfCurCell=index[i+1]-posOfCurCell;
+      INTERP_KERNEL::NormalizedCellType type=(INTERP_KERNEL::NormalizedCellType)conn[posOfCurCell];
+      int newLgth;
+      INTERP_KERNEL::NormalizedCellType newType=INTERP_KERNEL::CellSimplify::simplifyDegeneratedCell(type,conn+posOfCurCell+1,lgthOfCurCell-1,
+                                                                                                     conn+newPos+1,newLgth);
+      // Shall we delete the cell if it is completely degenerated:
+      bool delCell=INTERP_KERNEL::CellSimplify::isFlatCell(conn, newPos, newLgth, newType);
+      if (delCell)
+        {
+          nbDelCells++;
+          ret->pushBackSilent(i);
+        }
+      else   //if the cell is to be deleted, simply stay at the same place
+        {
+          conn[newPos]=newType;
+          newPos+=newLgth+1;
+        }
+      posOfCurCell=index[i+1];
+      index[i+1-nbDelCells]=newPos;
+    }
+  if(newPos!=initMeshLgth)
+    _nodal_connec->reAlloc(newPos);
+  const int nCellDel=ret->getNumberOfTuples();
+  if (nCellDel)
+    _nodal_connec_index->reAlloc(nbOfCells-nCellDel+1);
+  computeTypes();
+  return ret.retn();
+}
+
+/*!
+ * This method remove null 1D cells from \a this. A 1D cell is considered null if start node is equal to end node.
+ * Only connectivity is considered here.
+ */
+bool MEDCouplingUMesh::removeDegenerated1DCells()
+{
+  checkConnectivityFullyDefined();
+  if(getMeshDimension()!=1)
+    throw INTERP_KERNEL::Exception("MEDCouplingUMesh::removeDegenerated1DCells works on umeshes with meshdim equals to 1 !");
+  std::size_t nbCells(getNumberOfCells()),newSize(0),newSize2(0);
+  const int *conn(getNodalConnectivity()->begin()),*conni(getNodalConnectivityIndex()->begin());
+  {
+    for(std::size_t i=0;i<nbCells;i++)
+      {
+        INTERP_KERNEL::NormalizedCellType ct((INTERP_KERNEL::NormalizedCellType)conn[conni[i]]);
+        if(ct==INTERP_KERNEL::NORM_SEG2 || ct==INTERP_KERNEL::NORM_SEG3)
+          {
+            if(conn[conni[i]+1]!=conn[conni[i]+2])
+              {
+                newSize++;
+                newSize2+=conni[i+1]-conni[i];
+              }
+          }
+        else
+          {
+            std::ostringstream oss; oss << "MEDCouplingUMesh::removeDegenerated1DCells : cell #" << i << " in this is not of type SEG2/SEG3 !";
+            throw INTERP_KERNEL::Exception(oss.str());
+          }
+      }
+  }
+  if(newSize==nbCells)//no cells has been removed -> do nothing
+    return false;
+  MCAuto<DataArrayInt> newConn(DataArrayInt::New()),newConnI(DataArrayInt::New()); newConnI->alloc(newSize+1,1); newConn->alloc(newSize2,1);
+  int *newConnPtr(newConn->getPointer()),*newConnIPtr(newConnI->getPointer()); newConnIPtr[0]=0;
+  for(std::size_t i=0;i<nbCells;i++)
+    {
+      if(conn[conni[i]+1]!=conn[conni[i]+2])
+        {
+          newConnIPtr[1]=newConnIPtr[0]+conni[i+1]-conni[i];
+          newConnPtr=std::copy(conn+conni[i],conn+conni[i+1],newConnPtr);
+          newConnIPtr++;
+        }
+    }
+  setConnectivity(newConn,newConnI,true);
+  return true;
+}
+
 /*!
  * Finds incorrectly oriented cells of this 2D mesh in 3D space.
  * A cell is considered to be oriented correctly if an angle between its
  * normal vector and a given vector is less than \c PI / \c 2.
  *  \param [in] vec - 3 components of the vector specifying the correct orientation of
- *         cells. 
+ *         cells.
  *  \param [in] polyOnly - if \c true, only polygons are checked, else, all cells are
  *         checked.
  *  \param [in,out] cells - a vector returning ids of incorrectly oriented cells. It
@@ -4700,9 +4804,9 @@ void MEDCouplingUMesh::are2DCellsNotCorrectlyOriented(const double *vec, bool po
 /*!
  * Reverse connectivity of 2D cells whose orientation is not correct. A cell is
  * considered to be oriented correctly if an angle between its normal vector and a
- * given vector is less than \c PI / \c 2. 
+ * given vector is less than \c PI / \c 2.
  *  \param [in] vec - 3 components of the vector specifying the correct orientation of
- *         cells. 
+ *         cells.
  *  \param [in] polyOnly - if \c true, only polygons are checked, else, all cells are
  *         checked.
  *  \throw If \a this->getMeshDimension() != 2.
@@ -4810,7 +4914,7 @@ void MEDCouplingUMesh::arePolyhedronsNotCorrectlyOriented(std::vector<int>& cell
 
 /*!
  * Tries to fix connectivity of polyhedra, so that normal vector of all facets to point
- * out of the cell. 
+ * out of the cell.
  *  \throw If \a this->getMeshDimension() != 3.
  *  \throw If \a this->getSpaceDimension() != 3.
  *  \throw If the coordinates array is not set.
@@ -4852,7 +4956,7 @@ void MEDCouplingUMesh::orientCorrectlyPolyhedrons()
 }
 
 /*!
- * This method invert orientation of all cells in \a this. 
+ * This method invert orientation of all cells in \a this.
  * After calling this method the absolute value of measure of cells in \a this are the same than before calling.
  * This method only operates on the connectivity so coordinates are not touched at all.
  */
@@ -4879,7 +4983,7 @@ void MEDCouplingUMesh::invertOrientationOfAllCells()
  * pointing out of cell.
  *  \return DataArrayInt * - a new instance of DataArrayInt holding ids of fixed
  *         cells. The caller is to delete this array using decrRef() as it is no more
- *         needed. 
+ *         needed.
  *  \throw If \a this->getMeshDimension() != 3.
  *  \throw If \a this->getSpaceDimension() != 3.
  *  \throw If the coordinates array is not set.
@@ -4924,9 +5028,9 @@ DataArrayInt *MEDCouplingUMesh::findAndCorrectBadOriented3DExtrudedCells()
  * This method is a faster method to correct orientation of all 3D cells in \a this.
  * This method works only if \a this is a 3D mesh, that is to say a mesh with mesh dimension 3 and a space dimension 3.
  * This method makes the hypothesis that \a this a coherent that is to say MEDCouplingUMesh::checkConsistency should throw no exception.
- * 
+ *
  * \return a newly allocated int array with one components containing cell ids renumbered to fit the convention of MED (MED file and MEDCoupling)
- * \sa MEDCouplingUMesh::orientCorrectlyPolyhedrons, 
+ * \sa MEDCouplingUMesh::orientCorrectlyPolyhedrons,
  */
 DataArrayInt *MEDCouplingUMesh::findAndCorrectBadOriented3DCells()
 {
@@ -5012,13 +5116,13 @@ void MEDCouplingUMesh::getFastAveragePlaneOfThis(double *vec, double *pos) const
  * INTERP_KERNEL::NORM_TRI3, INTERP_KERNEL::NORM_QUAD4 and INTERP_KERNEL::NORM_TETRA4.
  * For a cell of other type an exception is thrown.
  * Space dimension of a 2D mesh can be either 2 or 3.
- * The Edge Ratio of a cell \f$t\f$ is: 
+ * The Edge Ratio of a cell \f$t\f$ is:
  *  \f$\frac{|t|_\infty}{|t|_0}\f$,
  *  where \f$|t|_\infty\f$ and \f$|t|_0\f$ respectively denote the greatest and
  *  the smallest edge lengths of \f$t\f$.
  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble on
  *          cells and one time, lying on \a this mesh. The caller is to delete this
- *          field using decrRef() as it is no more needed. 
+ *          field using decrRef() as it is no more needed.
  *  \throw If the coordinates array is not set.
  *  \throw If \a this mesh contains elements of dimension different from the mesh dimension.
  *  \throw If the connectivity data array has more than one component.
@@ -5090,7 +5194,7 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getEdgeRatioField() const
  * Space dimension of a 2D mesh can be either 2 or 3.
  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble on
  *          cells and one time, lying on \a this mesh. The caller is to delete this
- *          field using decrRef() as it is no more needed. 
+ *          field using decrRef() as it is no more needed.
  *  \throw If the coordinates array is not set.
  *  \throw If \a this mesh contains elements of dimension different from the mesh dimension.
  *  \throw If the connectivity data array has more than one component.
@@ -5171,7 +5275,7 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getAspectRatioField() const
  *  \f]
  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble on
  *          cells and one time, lying on \a this mesh. The caller is to delete this
- *          field using decrRef() as it is no more needed. 
+ *          field using decrRef() as it is no more needed.
  *  \throw If the coordinates array is not set.
  *  \throw If \a this mesh contains elements of dimension different from the mesh dimension.
  *  \throw If the connectivity data array has more than one component.
@@ -5238,7 +5342,7 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::getWarpField() const
  * For a cell of other type an exception is thrown.
  *  \return MEDCouplingFieldDouble * - a new instance of MEDCouplingFieldDouble on
  *          cells and one time, lying on \a this mesh. The caller is to delete this
- *          field using decrRef() as it is no more needed. 
+ *          field using decrRef() as it is no more needed.
  *  \throw If the coordinates array is not set.
  *  \throw If \a this mesh contains elements of dimension different from the mesh dimension.
  *  \throw If the connectivity data array has more than one component.
@@ -5320,11 +5424,11 @@ MEDCouplingFieldDouble *MEDCouplingUMesh::computeDiameterField() const
 
 /*!
  * This method aggregate the bbox of each cell and put it into bbox parameter (xmin,xmax,ymin,ymax,zmin,zmax).
- * 
+ *
  * \param [in] arcDetEps - a parameter specifying in case of 2D quadratic polygon cell the detection limit between linear and arc circle. (By default 1e-12)
  *                         For all other cases this input parameter is ignored.
  * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components.
- * 
+ *
  * \throw If \a this is not fully set (coordinates and connectivity).
  * \throw If a cell in \a this has no valid nodeId.
  * \sa MEDCouplingUMesh::getBoundingBoxForBBTreeFast, MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic
@@ -5356,9 +5460,9 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree(double arcDetEps) con
 /*!
  * This method aggregate the bbox of each cell only considering the nodes constituting each cell and put it into bbox parameter.
  * So meshes having quadratic cells the computed bounding boxes can be invalid !
- * 
+ *
  * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components.
- * 
+ *
  * \throw If \a this is not fully set (coordinates and connectivity).
  * \throw If a cell in \a this has no valid nodeId.
  */
@@ -5406,7 +5510,7 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTreeFast() const
  * useful for 2D meshes having quadratic cells
  * because for this type of cells getBoundingBoxForBBTreeFast method may return invalid bounding boxes (since it just considers
  * the two extremities of the arc of circle).
- * 
+ *
  * \param [in] arcDetEps - a parameter specifying in case of 2D quadratic polygon cell the detection limit between linear and arc circle. (By default 1e-12)
  * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components.
  * \throw If \a this is not fully defined.
@@ -5417,7 +5521,7 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTreeFast() const
 DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arcDetEps) const
 {
   checkFullyDefined();
-  INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(arcDetEps);
+  INTERP_KERNEL::QuadraticPlanarPrecision arcPrec(arcDetEps);
 
   int spaceDim(getSpaceDimension()),mDim(getMeshDimension()),nbOfCells(getNumberOfCells());
   if(spaceDim!=2 || mDim!=2)
@@ -5442,7 +5546,7 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arc
       else
         pol=INTERP_KERNEL::QuadraticPolygon::BuildArcCirclePolygon(nodes);
       INTERP_KERNEL::Bounds b; b.prepareForAggregation(); pol->fillBounds(b); delete pol;
-      bbox[0]=b.getXMin(); bbox[1]=b.getXMax(); bbox[2]=b.getYMin(); bbox[3]=b.getYMax(); 
+      bbox[0]=b.getXMin(); bbox[1]=b.getXMax(); bbox[2]=b.getYMin(); bbox[3]=b.getYMax();
     }
   return ret.retn();
 }
@@ -5452,7 +5556,7 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arc
  * useful for 2D meshes having quadratic cells
  * because for this type of cells getBoundingBoxForBBTreeFast method may return invalid bounding boxes (since it just considers
  * the two extremities of the arc of circle).
- * 
+ *
  * \param [in] arcDetEps - a parameter specifying in case of 2D quadratic polygon cell the detection limit between linear and arc circle. (By default 1e-12)
  * \return DataArrayDouble * - newly created object (to be managed by the caller) \a this number of cells tuples and 2*spacedim components.
  * \throw If \a this is not fully defined.
@@ -5466,7 +5570,7 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic(double arc
   int spaceDim(getSpaceDimension()),mDim(getMeshDimension()),nbOfCells(getNumberOfCells());
   if(spaceDim!=2 || mDim!=1)
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic : This method should be applied on mesh with mesh dimension equal to 1 and space dimension also equal to 2!");
-  INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(arcDetEps);
+  INTERP_KERNEL::QuadraticPlanarPrecision arcPrec(arcDetEps);
   MCAuto<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim);
   double *bbox(ret->getPointer());
   const double *coords(_coords->begin());
@@ -5523,9 +5627,9 @@ namespace MEDCouplingImpl
  * This method expects that \a this is sorted by types. If not an exception will be thrown.
  * This method returns in the same format as code (see MEDCouplingUMesh::checkTypeConsistencyAndContig or MEDCouplingUMesh::splitProfilePerType) how
  * \a this is composed in cell types.
- * The returned array is of size 3*n where n is the number of different types present in \a this. 
- * For every k in [0,n] ret[3*k+2]==-1 because it has no sense here. 
- * This parameter is kept only for compatibility with other methode listed above.
+ * The returned array is of size 3*n where n is the number of different types present in \a this.
+ * For every k in [0,n] ret[3*k+2]==-1 because it has no sense here.
+ * This parameter is kept only for compatibility with other method listed above.
  */
 std::vector<int> MEDCouplingUMesh::getDistributionOfTypes() const
 {
@@ -5567,10 +5671,10 @@ std::vector<int> MEDCouplingUMesh::getDistributionOfTypes() const
  * If it exists k so that 3*k geometric type is not in geometric types of this an exception will be thrown.
  * If it exists k so that 3*k geometric type exists but the number of consecutive cell types does not match,
  * an exception is thrown too.
- * 
+ *
  * If all geometric types in \a code are exactly those in \a this null pointer is returned.
- * If it exists a geometric type in \a this \b not in \a code \b no exception is thrown 
- * and a DataArrayInt instance is returned that the user has the responsability to deallocate.
+ * If it exists a geometric type in \a this \b not in \a code \b no exception is thrown
+ * and a DataArrayInt instance is returned that the user has the responsibility to deallocate.
  */
 DataArrayInt *MEDCouplingUMesh::checkTypeConsistencyAndContig(const std::vector<int>& code, const std::vector<const DataArrayInt *>& idsPerType) const
 {
@@ -5664,7 +5768,7 @@ DataArrayInt *MEDCouplingUMesh::checkTypeConsistencyAndContig(const std::vector<
  * This method is the opposite of MEDCouplingUMesh::checkTypeConsistencyAndContig method. Given a list of cells in \a profile it returns a list of sub-profiles sorted by geo type.
  * The result is put in the array \a idsPerType. In the returned parameter \a code, foreach i \a code[3*i+2] refers (if different from -1) to a location into the \a idsPerType.
  * This method has 1 input \a profile and 3 outputs \a code \a idsInPflPerType and \a idsPerType.
- * 
+ *
  * \param [in] profile
  * \param [out] code is a vector of size 3*n where n is the number of different geometric type in \a this \b reduced to the profile \a profile. \a code has exactly the same semantic than in MEDCouplingUMesh::checkTypeConsistencyAndContig method.
  * \param [out] idsInPflPerType is a vector of size of different geometric type in the subpart defined by \a profile of \a this ( equal to \a code.size()/3). For each i,
@@ -5673,7 +5777,7 @@ DataArrayInt *MEDCouplingUMesh::checkTypeConsistencyAndContig(const std::vector<
  *              This vector can be empty in case of all geometric type cells are fully covered in ascending in the given input \a profile.
  * \throw if \a profile has not exactly one component. It throws too, if \a profile contains some values not in [0,getNumberOfCells()) or if \a this is not fully defined
  */
-void MEDCouplingUMesh::splitProfilePerType(const DataArrayInt *profile, std::vector<int>& code, std::vector<DataArrayInt *>& idsInPflPerType, std::vector<DataArrayInt *>& idsPerType) const
+void MEDCouplingUMesh::splitProfilePerType(const DataArrayInt *profile, std::vector<int>& code, std::vector<DataArrayInt *>& idsInPflPerType, std::vector<DataArrayInt *>& idsPerType, bool smartPflKiller) const
 {
   if(!profile)
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::splitProfilePerType : input profile is NULL !");
@@ -5715,7 +5819,7 @@ void MEDCouplingUMesh::splitProfilePerType(const DataArrayInt *profile, std::vec
       code[3*i]=(int)types[castId];
       code[3*i+1]=tmp3->getNumberOfTuples();
       MCAuto<DataArrayInt> tmp4=rankInsideCast->selectByTupleId(tmp3->begin(),tmp3->begin()+tmp3->getNumberOfTuples());
-      if(!tmp4->isIota(typeRangeVals[castId+1]-typeRangeVals[castId]))
+      if(!smartPflKiller || !tmp4->isIota(typeRangeVals[castId+1]-typeRangeVals[castId]))
         {
           tmp4->copyStringInfoFrom(*profile);
           idsPerType2.push_back(tmp4);
@@ -5828,7 +5932,7 @@ bool MEDCouplingUMesh::checkConsecutiveCellTypes() const
 /*!
  * This method is a specialization of MEDCouplingUMesh::checkConsecutiveCellTypesAndOrder method that is called here.
  * The geometric type order is specified by MED file.
- * 
+ *
  * \sa  MEDCouplingUMesh::checkConsecutiveCellTypesAndOrder
  */
 bool MEDCouplingUMesh::checkConsecutiveCellTypesForMEDFileFrmt() const
@@ -5881,7 +5985,7 @@ bool MEDCouplingUMesh::checkConsecutiveCellTypesAndOrder(const INTERP_KERNEL::No
 /*!
  * This method returns 2 newly allocated DataArrayInt instances. The first is an array of size 'this->getNumberOfCells()' with one component,
  * that tells for each cell the pos of its type in the array on type given in input parameter. The 2nd output parameter is an array with the same
- * number of tuples than input type array and with one component. This 2nd output array gives type by type the number of occurence of type in 'this'.
+ * number of tuples than input type array and with one component. This 2nd output array gives type by type the number of occurrence of type in 'this'.
  */
 DataArrayInt *MEDCouplingUMesh::getLevArrPerCellTypes(const INTERP_KERNEL::NormalizedCellType *orderBg, const INTERP_KERNEL::NormalizedCellType *orderEnd, DataArrayInt *&nbPerType) const
 {
@@ -5920,7 +6024,7 @@ DataArrayInt *MEDCouplingUMesh::getLevArrPerCellTypes(const INTERP_KERNEL::Norma
 /*!
  * This method behaves exactly as MEDCouplingUMesh::getRenumArrForConsecutiveCellTypesSpec but the order is those defined in MED file spec.
  *
- * \return a new object containing the old to new correspondance.
+ * \return a new object containing the old to new correspondence.
  *
  * \sa MEDCouplingUMesh::getRenumArrForConsecutiveCellTypesSpec, MEDCouplingUMesh::sortCellsInMEDFileFrmt.
  */
@@ -5930,7 +6034,7 @@ DataArrayInt *MEDCouplingUMesh::getRenumArrForMEDFileFrmt() const
 }
 
 /*!
- * This method is similar to method MEDCouplingUMesh::rearrange2ConsecutiveCellTypes except that the type order is specfied by [ \a orderBg , \a orderEnd ) (as MEDCouplingUMesh::checkConsecutiveCellTypesAndOrder method) and that this method is \b const and performs \b NO permutation in \a this.
+ * This method is similar to method MEDCouplingUMesh::rearrange2ConsecutiveCellTypes except that the type order is specified by [ \a orderBg , \a orderEnd ) (as MEDCouplingUMesh::checkConsecutiveCellTypesAndOrder method) and that this method is \b const and performs \b NO permutation in \a this.
  * This method returns an array of size getNumberOfCells() that gives a renumber array old2New that can be used as input of MEDCouplingMesh::renumberCells.
  * The mesh after this call to MEDCouplingMesh::renumberCells will pass the test of MEDCouplingUMesh::checkConsecutiveCellTypesAndOrder with the same inputs.
  * The returned array minimizes the permutations that is to say the order of cells inside same geometric type remains the same.
@@ -5949,7 +6053,7 @@ DataArrayInt *MEDCouplingUMesh::getRenumArrForConsecutiveCellTypesSpec(const INT
  * This method tries to minimizes the number of needed permutations. So, this method behaves not exactly as
  * MEDCouplingUMesh::sortCellsInMEDFileFrmt.
  *
- * \return the array giving the correspondance old to new.
+ * \return the array giving the correspondence old to new.
  */
 DataArrayInt *MEDCouplingUMesh::rearrange2ConsecutiveCellTypes()
 {
@@ -6122,7 +6226,7 @@ void MEDCouplingUMesh::convertNodalConnectivityToDynamicGeoTypeMesh(DataArrayInt
 /*!
  * This method takes in input a vector of MEDCouplingUMesh instances lying on the same coordinates with same mesh dimensions.
  * Each mesh in \b ms must be sorted by type with the same order (typically using MEDCouplingUMesh::sortCellsInMEDFileFrmt).
- * This method is particulary useful for MED file interaction. It allows to aggregate several meshes and keeping the type sorting
+ * This method is particularly useful for MED file interaction. It allows to aggregate several meshes and keeping the type sorting
  * and the track of the permutation by chunk of same geotype cells to retrieve it. The traditional formats old2new and new2old
  * are not used here to avoid the build of big permutation array.
  *
@@ -6308,7 +6412,7 @@ MEDCouplingMesh *MEDCouplingUMesh::mergeMyselfWith(const MEDCouplingMesh *other)
 /*!
  * Returns a new DataArrayDouble holding barycenters of all cells. The barycenter is
  * computed by averaging coordinates of cell nodes, so this method is not a right
- * choice for degnerated meshes (not well oriented, cells with measure close to zero).
+ * choice for degenerated meshes (not well oriented, cells with measure close to zero).
  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size \a
  *          this->getNumberOfCells() tuples per \a this->getSpaceDimension()
  *          components. The caller is to delete this array using decrRef() as it is
@@ -6339,11 +6443,11 @@ DataArrayDouble *MEDCouplingUMesh::computeCellCenterOfMass() const
 
 /*!
  * This method computes for each cell in \a this, the location of the iso barycenter of nodes constituting
- * the cell. Contrary to badly named MEDCouplingUMesh::computeCellCenterOfMass method that returns the center of inertia of the 
- * 
- * \return a newly allocated DataArrayDouble instance that the caller has to deal with. The returned 
+ * the cell. Contrary to badly named MEDCouplingUMesh::computeCellCenterOfMass method that returns the center of inertia of the
+ *
+ * \return a newly allocated DataArrayDouble instance that the caller has to deal with. The returned
  *          DataArrayDouble instance will have \c this->getNumberOfCells() tuples and \c this->getSpaceDimension() components.
- * 
+ *
  * \sa MEDCouplingUMesh::computeCellCenterOfMass
  * \throw If \a this is not fully defined (coordinates and connectivity)
  * \throw If there is presence in nodal connectivity in \a this of node ids not in [0, \c this->getNumberOfNodes() )
@@ -6414,14 +6518,14 @@ DataArrayDouble *MEDCouplingUMesh::computeIsoBarycenterOfNodesPerCell() const
 /*!
  * Returns a new DataArrayDouble holding barycenters of specified cells. The
  * barycenter is computed by averaging coordinates of cell nodes. The cells to treat
- * are specified via an array of cell ids. 
- *  \warning Validity of the specified cell ids is not checked! 
+ * are specified via an array of cell ids.
+ *  \warning Validity of the specified cell ids is not checked!
  *           Valid range is [ 0, \a this->getNumberOfCells() ).
  *  \param [in] begin - an array of cell ids of interest.
  *  \param [in] end - the end of \a begin, i.e. a pointer to its (last+1)-th element.
  *  \return DataArrayDouble * - a new instance of DataArrayDouble, of size ( \a
  *          end - \a begin ) tuples per \a this->getSpaceDimension() components. The
- *          caller is to delete this array using decrRef() as it is no more needed. 
+ *          caller is to delete this array using decrRef() as it is no more needed.
  *  \throw If the coordinates array is not set.
  *  \throw If the nodal connectivity of cells is not defined.
  *
@@ -6457,7 +6561,7 @@ DataArrayDouble *MEDCouplingUMesh::getPartBarycenterAndOwner(const int *begin, c
  * So this method expects that \a this has a spaceDimension equal to 3 and meshDimension equal to 2.
  * The computation of the plane equation is done using each time the 3 first nodes of 2D cells.
  * This method is useful to detect 2D cells in 3D space that are not coplanar.
- * 
+ *
  * \return DataArrayDouble * - a new instance of DataArrayDouble having 4 components and a number of tuples equal to number of cells in \a this.
  * \throw If spaceDim!=3 or meshDim!=2.
  * \throw If connectivity of \a this is invalid.
@@ -6530,7 +6634,7 @@ DataArrayDouble *MEDCouplingUMesh::computePlaneEquationOf3DFaces() const
 
 /*!
  * This method expects as input a DataArrayDouble non nul instance 'da' that should be allocated. If not an exception is thrown.
- * 
+ *
  */
 MEDCouplingUMesh *MEDCouplingUMesh::Build0DMeshFromCoords(DataArrayDouble *da)
 {
@@ -6785,7 +6889,7 @@ MEDCouplingUMesh *MEDCouplingUMesh::FuseUMeshesOnSameCoords(const std::vector<co
  * Makes all given meshes share the nodal connectivity array. The common connectivity
  * array is created by concatenating the connectivity arrays of all given meshes. All
  * the given meshes must be of the same space dimension but dimension of cells **can
- * differ**. This method is particulary useful in MEDLoader context to build a \ref
+ * differ**. This method is particularly useful in MEDLoader context to build a \ref
  * MEDCoupling::MEDFileUMesh "MEDFileUMesh" instance that expects that underlying
  * MEDCouplingUMesh'es of different dimension share the same nodal connectivity array.
  *  \param [in,out] meshes - a vector of meshes to update.
@@ -6839,9 +6943,9 @@ void MEDCouplingUMesh::PutUMeshesOnSameAggregatedCoords(const std::vector<MEDCou
 /*!
  * Merges nodes coincident with a given precision within all given meshes that share
  * the nodal connectivity array. The given meshes **can be of different** mesh
- * dimension. This method is particulary useful in MEDLoader context to build a \ref
+ * dimension. This method is particularly useful in MEDLoader context to build a \ref
  * MEDCoupling::MEDFileUMesh "MEDFileUMesh" instance that expects that underlying
- * MEDCouplingUMesh'es of different dimension share the same nodal connectivity array. 
+ * MEDCouplingUMesh'es of different dimension share the same nodal connectivity array.
  *  \param [in,out] meshes - a vector of meshes to update.
  *  \param [in] eps - the precision used to detect coincident nodes (infinite norm).
  *  \throw If any of \a meshes is NULL.
@@ -6885,33 +6989,30 @@ void MEDCouplingUMesh::MergeNodesOnUMeshesSharingSameCoords(const std::vector<ME
     {
       (*it)->renumberNodesInConn(o2n->begin());
       (*it)->setCoords(newCoords);
-    } 
+    }
 }
 
 
 /*!
- * This static operates only for coords in 3D. The polygon is specfied by its connectivity nodes in [ \a begin , \a end ).
+ * This static operates only for coords in 3D. The polygon is specified by its connectivity nodes in [ \a begin , \a end ).
  */
 bool MEDCouplingUMesh::IsPolygonWellOriented(bool isQuadratic, const double *vec, const int *begin, const int *end, const double *coords)
 {
   std::size_t i, ip1;
   double v[3]={0.,0.,0.};
   std::size_t sz=std::distance(begin,end);
-  if(isQuadratic)
-    sz/=2;
-  for(i=0;i<sz;i++)
-    {
-      v[0]+=coords[3*begin[i]+1]*coords[3*begin[(i+1)%sz]+2]-coords[3*begin[i]+2]*coords[3*begin[(i+1)%sz]+1];
-      v[1]+=coords[3*begin[i]+2]*coords[3*begin[(i+1)%sz]]-coords[3*begin[i]]*coords[3*begin[(i+1)%sz]+2];
-      v[2]+=coords[3*begin[i]]*coords[3*begin[(i+1)%sz]+1]-coords[3*begin[i]+1]*coords[3*begin[(i+1)%sz]];
-    }
-  double ret = vec[0]*v[0]+vec[1]*v[1]+vec[2]*v[2];
-
-  // Try using quadratic points if standard points are degenerated (for example a QPOLYG with two
-  // SEG3 forming a circle):
-  if (fabs(ret) < INTERP_KERNEL::DEFAULT_ABS_TOL && isQuadratic)
+  if(!isQuadratic)
+    for(i=0;i<sz;i++)
+      {
+        v[0]+=coords[3*begin[i]+1]*coords[3*begin[(i+1)%sz]+2]-coords[3*begin[i]+2]*coords[3*begin[(i+1)%sz]+1];
+        v[1]+=coords[3*begin[i]+2]*coords[3*begin[(i+1)%sz]]-coords[3*begin[i]]*coords[3*begin[(i+1)%sz]+2];
+        v[2]+=coords[3*begin[i]]*coords[3*begin[(i+1)%sz]+1]-coords[3*begin[i]+1]*coords[3*begin[(i+1)%sz]];
+      }
+  else
     {
-      v[0] = 0.0; v[1] = 0.0; v[2] = 0.0;
+      // Use all points if quadratic (taking only linear points might lead to issues if the linearized version of the
+      // polygon is not convex or self-intersecting ... see testCellOrientation4)
+      sz /= 2;
       for(std::size_t j=0;j<sz;j++)
         {
           if (j%2)  // current point i is quadratic, next point i+1 is standard
@@ -6928,13 +7029,13 @@ bool MEDCouplingUMesh::IsPolygonWellOriented(bool isQuadratic, const double *vec
           v[1]+=coords[3*begin[i]+2]*coords[3*begin[ip1]]-coords[3*begin[i]]*coords[3*begin[ip1]+2];
           v[2]+=coords[3*begin[i]]*coords[3*begin[ip1]+1]-coords[3*begin[i]+1]*coords[3*begin[ip1]];
         }
-      ret = vec[0]*v[0]+vec[1]*v[1]+vec[2]*v[2];
     }
+  double ret = vec[0]*v[0]+vec[1]*v[1]+vec[2]*v[2];
   return (ret>0.);
 }
 
 /*!
- * The polyhedron is specfied by its connectivity nodes in [ \a begin , \a end ).
+ * The polyhedron is specified by its connectivity nodes in [ \a begin , \a end ).
  */
 bool MEDCouplingUMesh::IsPolyhedronWellOriented(const int *begin, const int *end, const double *coords)
 {
@@ -6994,7 +7095,7 @@ bool MEDCouplingUMesh::IsTetra4WellOriented(const int *begin, const int *end, co
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::IsTetra4WellOriented : Tetra4 cell with not 4 nodes ! Call checkConsistency !");
   double vec0[3],vec1[3];
   const double *pt0=coords+3*begin[0],*pt1=coords+3*begin[1],*pt2=coords+3*begin[2],*pt3=coords+3*begin[3];
-  vec0[0]=pt1[0]-pt0[0]; vec0[1]=pt1[1]-pt0[1]; vec0[2]=pt1[2]-pt0[2]; vec1[0]=pt2[0]-pt0[0]; vec1[1]=pt2[1]-pt0[1]; vec1[2]=pt2[2]-pt0[2]; 
+  vec0[0]=pt1[0]-pt0[0]; vec0[1]=pt1[1]-pt0[1]; vec0[2]=pt1[2]-pt0[2]; vec1[0]=pt2[0]-pt0[0]; vec1[1]=pt2[1]-pt0[1]; vec1[2]=pt2[2]-pt0[2];
   return ((vec0[1]*vec1[2]-vec0[2]*vec1[1])*(pt3[0]-pt0[0])+(vec0[2]*vec1[0]-vec0[0]*vec1[2])*(pt3[1]-pt0[1])+(vec0[0]*vec1[1]-vec0[1]*vec1[0])*(pt3[2]-pt0[2]))<0;
 }
 
@@ -7010,7 +7111,7 @@ bool MEDCouplingUMesh::IsPyra5WellOriented(const int *begin, const int *end, con
 }
 
 /*!
- * This method performs a simplyfication of a single polyedron cell. To do that each face of cell whose connectivity is defined by [ \b begin , \b end ) 
+ * This method performs a simplyfication of a single polyedron cell. To do that each face of cell whose connectivity is defined by [ \b begin , \b end )
  * is compared with the others in order to find faces in the same plane (with approx of eps). If any, the cells are grouped together and projected to
  * a 2D space.
  *
@@ -7020,48 +7121,62 @@ bool MEDCouplingUMesh::IsPyra5WellOriented(const int *begin, const int *end, con
  * \param [in] end end of nodal connectivity of a single polyhedron cell (excluded)
  * \param [out] res the result is put at the end of the vector without any alteration of the data.
  */
-void MEDCouplingUMesh::SimplifyPolyhedronCell(double eps, const DataArrayDouble *coords, const int *begin, const int *end, DataArrayInt *res)
+void MEDCouplingUMesh::SimplifyPolyhedronCell(double eps, const DataArrayDouble *coords, int index, DataArrayInt *res, MEDCouplingUMesh *faces,
+                                              DataArrayInt *E_Fi, DataArrayInt *E_F, DataArrayInt *F_Ei, DataArrayInt *F_E)
 {
-  int nbFaces=std::count(begin+1,end,-1)+1;
+  int nbFaces = E_Fi->getIJ(index + 1, 0) - E_Fi->getIJ(index, 0);
   MCAuto<DataArrayDouble> v=DataArrayDouble::New(); v->alloc(nbFaces,3);
   double *vPtr=v->getPointer();
-  MCAuto<DataArrayDouble> p=DataArrayDouble::New(); p->alloc(nbFaces,1);
+  MCAuto<DataArrayDouble> p=DataArrayDouble::New(); p->alloc(nbFaces,2);
   double *pPtr=p->getPointer();
-  const int *stFaceConn=begin+1;
+  int *e_fi = E_Fi->getPointer(), *e_f = E_F->getPointer(), *f_ei = F_Ei->getPointer(), *f_e = F_E->getPointer();
+  const int *f_idx = faces->getNodalConnectivityIndex()->getPointer(), *f_cnn = faces->getNodalConnectivity()->getPointer();
   for(int i=0;i<nbFaces;i++,vPtr+=3,pPtr++)
     {
-      const int *endFaceConn=std::find(stFaceConn,end,-1);
-      ComputeVecAndPtOfFace(eps,coords->begin(),stFaceConn,endFaceConn,vPtr,pPtr);
-      stFaceConn=endFaceConn+1;
+      int face = e_f[e_fi[index] + i];
+      ComputeVecAndPtOfFace(eps, coords->begin(), f_cnn + f_idx[face] + 1, f_cnn + f_idx[face + 1], vPtr, pPtr);
+      // to differentiate faces going to different cells:
+      pPtr++, *pPtr = 0;
+      for (int j = f_ei[face]; j < f_ei[face + 1]; j++)
+        *pPtr += f_e[j];
     }
   pPtr=p->getPointer(); vPtr=v->getPointer();
   DataArrayInt *comm1=0,*commI1=0;
   v->findCommonTuples(eps,-1,comm1,commI1);
+  for (int i = 0; i < nbFaces; i++)
+    if (comm1->findIdFirstEqual(i) < 0)
+      {
+        comm1->pushBackSilent(i);
+        commI1->pushBackSilent(comm1->getNumberOfTuples());
+      }
   MCAuto<DataArrayInt> comm1Auto(comm1),commI1Auto(commI1);
   const int *comm1Ptr=comm1->begin();
   const int *commI1Ptr=commI1->begin();
   int nbOfGrps1=commI1Auto->getNumberOfTuples()-1;
   res->pushBackSilent((int)INTERP_KERNEL::NORM_POLYHED);
   //
-  MCAuto<MEDCouplingUMesh> mm=MEDCouplingUMesh::New("",3);
-  mm->setCoords(const_cast<DataArrayDouble *>(coords)); mm->allocateCells(1); mm->insertNextCell(INTERP_KERNEL::NORM_POLYHED,(int)std::distance(begin+1,end),begin+1);
-  mm->finishInsertingCells();
-  //
   for(int i=0;i<nbOfGrps1;i++)
     {
       int vecId=comm1Ptr[commI1Ptr[i]];
       MCAuto<DataArrayDouble> tmpgrp2=p->selectByTupleId(comm1Ptr+commI1Ptr[i],comm1Ptr+commI1Ptr[i+1]);
       DataArrayInt *comm2=0,*commI2=0;
       tmpgrp2->findCommonTuples(eps,-1,comm2,commI2);
+      for (int j = 0; j < commI1Ptr[i+1] - commI1Ptr[i]; j++)
+        if (comm2->findIdFirstEqual(j) < 0)
+          {
+            comm2->pushBackSilent(j);
+            commI2->pushBackSilent(comm2->getNumberOfTuples());
+          }
       MCAuto<DataArrayInt> comm2Auto(comm2),commI2Auto(commI2);
       const int *comm2Ptr=comm2->begin();
       const int *commI2Ptr=commI2->begin();
       int nbOfGrps2=commI2Auto->getNumberOfTuples()-1;
       for(int j=0;j<nbOfGrps2;j++)
         {
-          if(commI2Ptr[j+1]-commI2Ptr[j]<=1)
+          if(commI2Ptr[j+1] == commI2Ptr[j] + 1)
             {
-              res->insertAtTheEnd(begin,end);
+              int face = e_f[e_fi[index] + comm1Ptr[commI1Ptr[i] + comm2Ptr[commI2Ptr[j]]]]; //hmmm
+              res->insertAtTheEnd(f_cnn + f_idx[face] + 1, f_cnn + f_idx[face + 1]);
               res->pushBackSilent(-1);
             }
           else
@@ -7069,13 +7184,12 @@ void MEDCouplingUMesh::SimplifyPolyhedronCell(double eps, const DataArrayDouble
               int pointId=comm1Ptr[commI1Ptr[i]+comm2Ptr[commI2Ptr[j]]];
               MCAuto<DataArrayInt> ids2=comm2->selectByTupleIdSafeSlice(commI2Ptr[j],commI2Ptr[j+1],1);
               ids2->transformWithIndArr(comm1Ptr+commI1Ptr[i],comm1Ptr+commI1Ptr[i+1]);
-              DataArrayInt *tmp0=DataArrayInt::New(),*tmp1=DataArrayInt::New(),*tmp2=DataArrayInt::New(),*tmp3=DataArrayInt::New();
-              MCAuto<MEDCouplingUMesh> mm2=mm->buildDescendingConnectivity(tmp0,tmp1,tmp2,tmp3); tmp0->decrRef(); tmp1->decrRef(); tmp2->decrRef(); tmp3->decrRef();
-              MCAuto<MEDCouplingUMesh> mm3=static_cast<MEDCouplingUMesh *>(mm2->buildPartOfMySelf(ids2->begin(),ids2->end(),true));
+              ids2->transformWithIndArr(e_f + e_fi[index], e_f + e_fi[index + 1]);
+              MCAuto<MEDCouplingUMesh> mm3=static_cast<MEDCouplingUMesh *>(faces->buildPartOfMySelf(ids2->begin(),ids2->end(),true));
               MCAuto<DataArrayInt> idsNodeTmp=mm3->zipCoordsTraducer();
               MCAuto<DataArrayInt> idsNode=idsNodeTmp->invertArrayO2N2N2O(mm3->getNumberOfNodes());
               const int *idsNodePtr=idsNode->begin();
-              double center[3]; center[0]=pPtr[pointId]*vPtr[3*vecId]; center[1]=pPtr[pointId]*vPtr[3*vecId+1]; center[2]=pPtr[pointId]*vPtr[3*vecId+2];
+              double center[3]; center[0]=pPtr[2*pointId]*vPtr[3*vecId]; center[1]=pPtr[2*pointId]*vPtr[3*vecId+1]; center[2]=pPtr[2*pointId]*vPtr[3*vecId+2];
               double vec[3]; vec[0]=vPtr[3*vecId+1]; vec[1]=-vPtr[3*vecId]; vec[2]=0.;
               double norm=vec[0]*vec[0]+vec[1]*vec[1]+vec[2]*vec[2];
               if(std::abs(norm)>eps)
@@ -7104,7 +7218,7 @@ void MEDCouplingUMesh::SimplifyPolyhedronCell(double eps, const DataArrayDouble
 /*!
  * This method computes the normalized vector of the plane and the pos of the point belonging to the plane and the line defined by the vector going
  * through origin. The plane is defined by its nodal connectivity [ \b begin, \b end ).
- * 
+ *
  * \param [in] eps below that value the dot product of 2 vectors is considered as colinears
  * \param [in] coords coordinates expected to have 3 components.
  * \param [in] begin start of the nodal connectivity of the face.
@@ -7188,7 +7302,7 @@ void MEDCouplingUMesh::TryToCorrectPolyhedronOrientation(int *begin, int *end, c
                   if(b1 || b2) { b=b2; isPerm[i]=true; smthChanged++; break; }
                 }
               if(isPerm[i])
-                { 
+                {
                   if(!b)
                     std::reverse(bgFace+1,endFace);
                   for(std::size_t j=0;j<nbOfEdgesInFace;j++)
@@ -7233,7 +7347,7 @@ void MEDCouplingUMesh::TryToCorrectPolyhedronOrientation(int *begin, int *end, c
 /*!
  * This method makes the assumption spacedimension == meshdimension == 2.
  * This method works only for linear cells.
- * 
+ *
  * \return a newly allocated array containing the connectivity of a polygon type enum included (NORM_POLYGON in pos#0)
  */
 DataArrayInt *MEDCouplingUMesh::buildUnionOf2DMesh() const
@@ -7257,7 +7371,7 @@ DataArrayInt *MEDCouplingUMesh::buildUnionOf2DMesh() const
 /*!
  * This method makes the assumption spacedimension == meshdimension == 3.
  * This method works only for linear cells.
- * 
+ *
  * \return a newly allocated array containing the connectivity of a polygon type enum included (NORM_POLYHED in pos#0)
  */
 DataArrayInt *MEDCouplingUMesh::buildUnionOf3DMesh() const
@@ -7568,7 +7682,7 @@ DataArrayInt *MEDCouplingUMesh::orderConsecutiveCells1D() const
  * avoid to have a non conform mesh.
  *
  * \return int - the number of new nodes created (in most of cases 0).
- * 
+ *
  * \throw If \a this is not coherent.
  * \throw If \a this has not spaceDim equal to 2.
  * \throw If \a this has not meshDim equal to 2.
@@ -7599,7 +7713,7 @@ int MEDCouplingUMesh::split2DCells(const DataArrayInt *desc, const DataArrayInt
  * the geometric cell type set to INTERP_KERNEL::NORM_POLYGON.
  * This method excepts that \b coords parameter is expected to be in dimension 2. [ \b nodalConnBg , \b nodalConnEnd ) is the nodal connectivity of the input
  * cell (geometric cell type included at the position 0). If the meshdimension of the input cell is not equal to 2 an INTERP_KERNEL::Exception will be thrown.
- * 
+ *
  * \return false if the input connectivity represents already the convex hull, true if the input cell needs to be reordered.
  */
 bool MEDCouplingUMesh::BuildConvexEnvelopOf2DCellJarvis(const double *coords, const int *nodalConnBg, const int *nodalConnEnd, DataArrayInt *nodalConnecOut)
@@ -7693,7 +7807,7 @@ bool MEDCouplingUMesh::BuildConvexEnvelopOf2DCellJarvis(const double *coords, co
 /*!
  * This method works on an input pair (\b arr, \b arrIndx) where \b arr indexes is in \b arrIndx.
  * This method will not impact the size of inout parameter \b arrIndx but the size of \b arr will be modified in case of suppression.
- * 
+ *
  * \param [in] idsToRemoveBg begin of set of ids to remove in \b arr (included)
  * \param [in] idsToRemoveEnd end of set of ids to remove in \b arr (excluded)
  * \param [in,out] arr array in which the remove operation will be done.
@@ -7885,7 +7999,7 @@ void MEDCouplingUMesh::ExtractFromIndexedArraysSlice(int idsOfSelectStart, int i
  * This method works on an input pair (\b arrIn, \b arrIndxIn) where \b arrIn indexes is in \b arrIndxIn.
  * This method builds an output pair (\b arrOut,\b arrIndexOut) that is a copy from \b arrIn for all cell ids \b not \b in [ \b idsOfSelectBg , \b idsOfSelectEnd ) and for
  * cellIds \b in [ \b idsOfSelectBg , \b idsOfSelectEnd ) a copy coming from the corresponding values in input pair (\b srcArr, \b srcArrIndex).
- * This method is an generalization of MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx that performs the same thing but by without building explicitely a result output arrays.
+ * This method is an generalization of MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx that performs the same thing but by without building explicitly a result output arrays.
  *
  * \param [in] idsOfSelectBg begin of set of ids of the input extraction (included)
  * \param [in] idsOfSelectEnd end of set of ids of the input extraction (excluded)
@@ -7895,7 +8009,7 @@ void MEDCouplingUMesh::ExtractFromIndexedArraysSlice(int idsOfSelectStart, int i
  * \param [in] srcArrIndex index array of \b srcArr
  * \param [out] arrOut the resulting array
  * \param [out] arrIndexOut the index array of the resulting array \b arrOut
- * 
+ *
  * \sa MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx
  */
 void MEDCouplingUMesh::SetPartOfIndexedArrays(const int *idsOfSelectBg, const int *idsOfSelectEnd, const DataArrayInt *arrIn, const DataArrayInt *arrIndxIn,
@@ -7951,7 +8065,7 @@ void MEDCouplingUMesh::SetPartOfIndexedArrays(const int *idsOfSelectBg, const in
 
 /*!
  * This method works on an input pair (\b arrIn, \b arrIndxIn) where \b arrIn indexes is in \b arrIndxIn.
- * This method is an specialization of MEDCouplingUMesh::SetPartOfIndexedArrays in the case of assignement do not modify the index in \b arrIndxIn.
+ * This method is an specialization of MEDCouplingUMesh::SetPartOfIndexedArrays in the case of assignment do not modify the index in \b arrIndxIn.
  *
  * \param [in] idsOfSelectBg begin of set of ids of the input extraction (included)
  * \param [in] idsOfSelectEnd end of set of ids of the input extraction (excluded)
@@ -7959,7 +8073,7 @@ void MEDCouplingUMesh::SetPartOfIndexedArrays(const int *idsOfSelectBg, const in
  * \param [in] arrIndxIn is the input index array allowing to walk into \b arrIn
  * \param [in] srcArr input array that will be used as source of copy for ids in [ \b idsOfSelectBg , \b idsOfSelectEnd )
  * \param [in] srcArrIndex index array of \b srcArr
- * 
+ *
  * \sa MEDCouplingUMesh::SetPartOfIndexedArrays
  */
 void MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx(const int *idsOfSelectBg, const int *idsOfSelectEnd, DataArrayInt *arrInOut, const DataArrayInt *arrIndxIn,
@@ -7999,7 +8113,7 @@ void MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx(const int *idsOfSelectBg, c
  * Then it is repeated recursively until either all ids are fetched or no more ids are reachable step by step.
  * A negative value in \b arrIn means that it is ignored.
  * This method is useful to see if a mesh is contiguous regarding its connectivity. If it is not the case the size of returned array is different from arrIndxIn->getNumberOfTuples()-1.
- * 
+ *
  * \param [in] arrIn arr origin array from which the extraction will be done.
  * \param [in] arrIndxIn is the input index array allowing to walk into \b arrIn
  * \return a newly allocated DataArray that stores all ids fetched by the gradually spread process.
@@ -8048,7 +8162,7 @@ DataArrayInt *MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeed(const int *se
  * This method works on an input pair (\b arrIn, \b arrIndxIn) where \b arrIn indexes is in \b arrIndxIn.
  * This method builds an output pair (\b arrOut,\b arrIndexOut) that is a copy from \b arrIn for all cell ids \b not \b in [ \b idsOfSelectBg , \b idsOfSelectEnd ) and for
  * cellIds \b in [\b idsOfSelectBg, \b idsOfSelectEnd) a copy coming from the corresponding values in input pair (\b srcArr, \b srcArrIndex).
- * This method is an generalization of MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx that performs the same thing but by without building explicitely a result output arrays.
+ * This method is an generalization of MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx that performs the same thing but by without building explicitly a result output arrays.
  *
  * \param [in] start begin of set of ids of the input extraction (included)
  * \param [in] end end of set of ids of the input extraction (excluded)
@@ -8059,7 +8173,7 @@ DataArrayInt *MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeed(const int *se
  * \param [in] srcArrIndex index array of \b srcArr
  * \param [out] arrOut the resulting array
  * \param [out] arrIndexOut the index array of the resulting array \b arrOut
- * 
+ *
  * \sa MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx MEDCouplingUMesh::SetPartOfIndexedArrays
  */
 void MEDCouplingUMesh::SetPartOfIndexedArraysSlice(int start, int end, int step, const DataArrayInt *arrIn, const DataArrayInt *arrIndxIn,
@@ -8113,7 +8227,7 @@ void MEDCouplingUMesh::SetPartOfIndexedArraysSlice(int start, int end, int step,
 
 /*!
  * This method works on an input pair (\b arrIn, \b arrIndxIn) where \b arrIn indexes is in \b arrIndxIn.
- * This method is an specialization of MEDCouplingUMesh::SetPartOfIndexedArrays in the case of assignement do not modify the index in \b arrIndxIn.
+ * This method is an specialization of MEDCouplingUMesh::SetPartOfIndexedArrays in the case of assignment do not modify the index in \b arrIndxIn.
  *
  * \param [in] start begin of set of ids of the input extraction (included)
  * \param [in] end end of set of ids of the input extraction (excluded)
@@ -8122,7 +8236,7 @@ void MEDCouplingUMesh::SetPartOfIndexedArraysSlice(int start, int end, int step,
  * \param [in] arrIndxIn is the input index array allowing to walk into \b arrIn
  * \param [in] srcArr input array that will be used as source of copy for ids in [\b idsOfSelectBg, \b idsOfSelectEnd)
  * \param [in] srcArrIndex index array of \b srcArr
- * 
+ *
  * \sa MEDCouplingUMesh::SetPartOfIndexedArraysSlice MEDCouplingUMesh::SetPartOfIndexedArraysSameIdx
  */
 void MEDCouplingUMesh::SetPartOfIndexedArraysSameIdxSlice(int start, int end, int step, DataArrayInt *arrInOut, const DataArrayInt *arrIndxIn,
@@ -8163,7 +8277,7 @@ void MEDCouplingUMesh::SetPartOfIndexedArraysSameIdxSlice(int start, int end, in
  * The returned mesh contains as poly cells as number of contiguous zone (regarding connectivity).
  * A spread contiguous zone is built using poly cells (polyhedra in 3D, polygons in 2D and polyline in 1D).
  * The sum of measure field of returned mesh is equal to the sum of measure field of this.
- * 
+ *
  * \return a newly allocated mesh lying on the same coords than \b this with same meshdimension than \b this.
  */
 MEDCouplingUMesh *MEDCouplingUMesh::buildSpreadZonesWithPoly() const
@@ -8272,20 +8386,20 @@ DataArrayInt *MEDCouplingUMesh::ComputeRangesFromTypeDistribution(const std::vec
  * This method will split **all** 3D cells in \a this into INTERP_KERNEL::NORM_TETRA4 cells and put them in the returned mesh.
  * It leads to an increase to number of cells.
  * This method contrary to MEDCouplingUMesh::simplexize can append coordinates in \a this to perform its work.
- * The \a nbOfAdditionalPoints returned value informs about it. If > 0, the coordinates array in returned mesh will have \a nbOfAdditionalPoints 
+ * The \a nbOfAdditionalPoints returned value informs about it. If > 0, the coordinates array in returned mesh will have \a nbOfAdditionalPoints
  * more tuples (nodes) than in \a this. Anyway, all the nodes in \a this (with the same order) will be in the returned mesh.
  *
  * \param [in] policy - the policy of splitting that must be in (PLANAR_FACE_5, PLANAR_FACE_6, GENERAL_24, GENERAL_48). The policy will be used only for INTERP_KERNEL::NORM_HEXA8 cells.
  *                      For all other cells, the splitting policy will be ignored. See INTERP_KERNEL::SplittingPolicy for the images.
- * \param [out] nbOfAdditionalPoints - number of nodes added to \c this->_coords. If > 0 a new coordinates object will be constructed result of the aggregation of the old one and the new points added. 
+ * \param [out] nbOfAdditionalPoints - number of nodes added to \c this->_coords. If > 0 a new coordinates object will be constructed result of the aggregation of the old one and the new points added.
  * \param [out] n2oCells - A new instance of DataArrayInt holding, for each new cell,
  *          an id of old cell producing it. The caller is to delete this array using
  *         decrRef() as it is no more needed.
  * \return MEDCoupling1SGTUMesh * - the mesh containing only INTERP_KERNEL::NORM_TETRA4 cells.
  *
- * \warning This method operates on each cells in this independantly ! So it can leads to non conform mesh in returned value ! If you expect to have a conform mesh in output
+ * \warning This method operates on each cells in this independently ! So it can leads to non conform mesh in returned value ! If you expect to have a conform mesh in output
  * the policy PLANAR_FACE_6 should be used on a mesh sorted with MEDCoupling1SGTUMesh::sortHexa8EachOther.
- * 
+ *
  * \throw If \a this is not a 3D mesh (spaceDim==3 and meshDim==3).
  * \throw If \a this is not fully constituted with linear 3D cells.
  * \sa MEDCouplingUMesh::simplexize, MEDCoupling1SGTUMesh::sortHexa8EachOther