Salome HOME
scotch6.0.4 needs pthread... Quick and dirty solution to be improved
[tools/medcoupling.git] / src / MEDCoupling / MEDCouplingUMesh.cxx
index 907b00d6f09955cc249b0bc5dba8995480abc691..ee74930096d04567a3a9f28742b0e043d680a8e2 100644 (file)
@@ -16,7 +16,7 @@
 //
 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
-// Author : Anthony Geay (CEA/DEN)
+// Author : Anthony Geay (EDF R&D)
 
 #include "MEDCouplingUMesh.hxx"
 #include "MEDCouplingCMesh.hxx"
@@ -41,6 +41,7 @@
 #include "InterpKernelGeo2DEdgeLin.hxx"
 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
+#include "OrientationInverter.hxx"
 #include "MEDCouplingUMesh_internal.hxx"
 
 #include <sstream>
@@ -55,8 +56,8 @@ using namespace MEDCoupling;
 double MEDCouplingUMesh::EPS_FOR_POLYH_ORIENTATION=1.e-14;
 
 /// @cond INTERNAL
-const INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::MEDMEM_ORDER[N_MEDMEM_ORDER] = { INTERP_KERNEL::NORM_POINT1, INTERP_KERNEL::NORM_SEG2, INTERP_KERNEL::NORM_SEG3, INTERP_KERNEL::NORM_SEG4, INTERP_KERNEL::NORM_POLYL, INTERP_KERNEL::NORM_TRI3, INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6, INTERP_KERNEL::NORM_TRI7, INTERP_KERNEL::NORM_QUAD8, INTERP_KERNEL::NORM_QUAD9, INTERP_KERNEL::NORM_POLYGON, INTERP_KERNEL::NORM_QPOLYG, INTERP_KERNEL::NORM_TETRA4, INTERP_KERNEL::NORM_PYRA5, INTERP_KERNEL::NORM_PENTA6, INTERP_KERNEL::NORM_HEXA8, INTERP_KERNEL::NORM_HEXGP12, INTERP_KERNEL::NORM_TETRA10, INTERP_KERNEL::NORM_PYRA13, INTERP_KERNEL::NORM_PENTA15, INTERP_KERNEL::NORM_HEXA20, INTERP_KERNEL::NORM_HEXA27, INTERP_KERNEL::NORM_POLYHED };
-const int MEDCouplingUMesh::MEDCOUPLING2VTKTYPETRADUCER[INTERP_KERNEL::NORM_MAXTYPE+1]={1,3,21,5,9,7,22,34,23,28,-1,-1,-1,-1,10,14,13,-1,12,-1,24,-1,16,27,-1,26,-1,29,-1,-1,25,42,36,4};
+const INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::MEDMEM_ORDER[N_MEDMEM_ORDER] = { INTERP_KERNEL::NORM_POINT1, INTERP_KERNEL::NORM_SEG2, INTERP_KERNEL::NORM_SEG3, INTERP_KERNEL::NORM_SEG4, INTERP_KERNEL::NORM_POLYL, INTERP_KERNEL::NORM_TRI3, INTERP_KERNEL::NORM_QUAD4, INTERP_KERNEL::NORM_TRI6, INTERP_KERNEL::NORM_TRI7, INTERP_KERNEL::NORM_QUAD8, INTERP_KERNEL::NORM_QUAD9, INTERP_KERNEL::NORM_POLYGON, INTERP_KERNEL::NORM_QPOLYG, INTERP_KERNEL::NORM_TETRA4, INTERP_KERNEL::NORM_PYRA5, INTERP_KERNEL::NORM_PENTA6, INTERP_KERNEL::NORM_HEXA8, INTERP_KERNEL::NORM_HEXGP12, INTERP_KERNEL::NORM_TETRA10, INTERP_KERNEL::NORM_PYRA13, INTERP_KERNEL::NORM_PENTA15, INTERP_KERNEL::NORM_PENTA18, INTERP_KERNEL::NORM_HEXA20, INTERP_KERNEL::NORM_HEXA27, INTERP_KERNEL::NORM_POLYHED };
+const int MEDCouplingUMesh::MEDCOUPLING2VTKTYPETRADUCER[INTERP_KERNEL::NORM_MAXTYPE+1]={1,3,21,5,9,7,22,34,23,28,-1,-1,-1,-1,10,14,13,-1,12,-1,24,-1,16,27,-1,26,-1,29,32,-1,25,42,36,4};
 /// @endcond
 
 MEDCouplingUMesh *MEDCouplingUMesh::New()
@@ -300,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.
  *
@@ -589,18 +591,15 @@ void MEDCouplingUMesh::checkFastEquivalWith(const MEDCouplingMesh *other, double
 void MEDCouplingUMesh::getReverseNodalConnectivity(DataArrayInt *revNodal, DataArrayInt *revNodalIndx) const
 {
   checkFullyDefined();
-  int nbOfNodes=getNumberOfNodes();
+  int nbOfNodes(getNumberOfNodes());
   int *revNodalIndxPtr=(int *)malloc((nbOfNodes+1)*sizeof(int));
   revNodalIndx->useArray(revNodalIndxPtr,true,C_DEALLOC,nbOfNodes+1,1);
   std::fill(revNodalIndxPtr,revNodalIndxPtr+nbOfNodes+1,0);
-  const int *conn=_nodal_connec->getConstPointer();
-  const int *connIndex=_nodal_connec_index->getConstPointer();
-  int nbOfCells=getNumberOfCells();
-  int nbOfEltsInRevNodal=0;
+  const int *conn(_nodal_connec->begin()),*connIndex(_nodal_connec_index->begin());
+  int nbOfCells(getNumberOfCells()),nbOfEltsInRevNodal(0);
   for(int eltId=0;eltId<nbOfCells;eltId++)
     {
-      const int *strtNdlConnOfCurCell=conn+connIndex[eltId]+1;
-      const int *endNdlConnOfCurCell=conn+connIndex[eltId+1];
+      const int *strtNdlConnOfCurCell(conn+connIndex[eltId]+1),*endNdlConnOfCurCell(conn+connIndex[eltId+1]);
       for(const int *iter=strtNdlConnOfCurCell;iter!=endNdlConnOfCurCell;iter++)
         if(*iter>=0)//for polyhedrons
           {
@@ -838,10 +837,10 @@ void MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(const DataArrayInt *desc, cons
 {
   if(!desc || !descIndx || !revDesc || !revDescIndx)
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::ComputeNeighborsOfCellsAdv some input array is empty !");
-  const int *descPtr=desc->getConstPointer();
-  const int *descIPtr=descIndx->getConstPointer();
-  const int *revDescPtr=revDesc->getConstPointer();
-  const int *revDescIPtr=revDescIndx->getConstPointer();
+  const int *descPtr=desc->begin();
+  const int *descIPtr=descIndx->begin();
+  const int *revDescPtr=revDesc->begin();
+  const int *revDescIPtr=revDescIndx->begin();
   //
   int nbCells=descIndx->getNumberOfTuples()-1;
   MCAuto<DataArrayInt> out0=DataArrayInt::New();
@@ -863,6 +862,35 @@ void MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(const DataArrayInt *desc, cons
   neighborsIndx=out1.retn();
 }
 
+/*!
+ * Explodes \a this into edges whatever its dimension.
+ */
+MCAuto<MEDCouplingUMesh> MEDCouplingUMesh::explodeIntoEdges(MCAuto<DataArrayInt>& desc, MCAuto<DataArrayInt>& descIndex, MCAuto<DataArrayInt>& revDesc, MCAuto<DataArrayInt>& revDescIndx) const
+{
+  checkFullyDefined();
+  int mdim(getMeshDimension());
+  desc=DataArrayInt::New(); descIndex=DataArrayInt::New(); revDesc=DataArrayInt::New(); revDescIndx=DataArrayInt::New();
+  MCAuto<MEDCouplingUMesh> mesh1D;
+  switch(mdim)
+  {
+    case 3:
+      {
+        mesh1D=explode3DMeshTo1D(desc,descIndex,revDesc,revDescIndx);
+        break;
+      }
+    case 2:
+      {
+        mesh1D=buildDescendingConnectivity(desc,descIndex,revDesc,revDescIndx);
+        break;
+      }
+    default:
+      {
+        throw INTERP_KERNEL::Exception("MEDCouplingUMesh::computeNeighborsOfNodes : Mesh dimension supported are [3,2] !");
+      }
+  }
+  return mesh1D;
+}
+
 /*!
  * \b WARNING this method do the assumption that connectivity lies on the coordinates set.
  * For speed reasons no check of this will be done. This method calls
@@ -877,13 +905,15 @@ void MEDCouplingUMesh::ComputeNeighborsOfCellsAdv(const DataArrayInt *desc, cons
  * 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
 {
   checkFullyDefined();
   int mdim(getMeshDimension()),nbNodes(getNumberOfNodes());
   MCAuto<DataArrayInt> desc(DataArrayInt::New()),descIndx(DataArrayInt::New()),revDesc(DataArrayInt::New()),revDescIndx(DataArrayInt::New());
-  MCAuto<MEDCouplingUMesh> mesh1D;
+  MCConstAuto<MEDCouplingUMesh> mesh1D;
   switch(mdim)
   {
     case 3:
@@ -898,8 +928,7 @@ void MEDCouplingUMesh::computeNeighborsOfNodes(DataArrayInt *&neighbors, DataArr
       }
     case 1:
       {
-        mesh1D=const_cast<MEDCouplingUMesh *>(this);
-        mesh1D->incrRef();
+        mesh1D.takeRef(this);
         break;
       }
     default:
@@ -922,6 +951,50 @@ void MEDCouplingUMesh::computeNeighborsOfNodes(DataArrayInt *&neighbors, DataArr
   neighborsIdx=descIndx.retn();
 }
 
+/*!
+ * 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
+{
+  checkFullyDefined();
+  int nbOfNodes(getNumberOfNodes());
+  const int *conn(_nodal_connec->begin()),*connIndex(_nodal_connec_index->begin());
+  int nbOfCells(getNumberOfCells());
+  std::vector< std::set<int> > st0(nbOfNodes);
+  for(int eltId=0;eltId<nbOfCells;eltId++)
+    {
+      const int *strtNdlConnOfCurCell(conn+connIndex[eltId]+1),*endNdlConnOfCurCell(conn+connIndex[eltId+1]);
+      std::set<int> s(strtNdlConnOfCurCell,endNdlConnOfCurCell); s.erase(-1); //for polyhedrons
+      for(std::set<int>::const_iterator iter2=s.begin();iter2!=s.end();iter2++)
+        st0[*iter2].insert(s.begin(),s.end());
+    }
+  neighborsIdx=DataArrayInt::New(); neighborsIdx->alloc(nbOfNodes+1,1); neighborsIdx->setIJ(0,0,0);
+  {
+    int *neighIdx(neighborsIdx->getPointer());
+    for(std::vector< std::set<int> >::const_iterator it=st0.begin();it!=st0.end();it++,neighIdx++)
+      {
+        if ((*it).empty())
+          neighIdx[1]=neighIdx[0];
+        else
+          neighIdx[1]=neighIdx[0]+(*it).size()-1;
+      }
+  }
+  neighbors=DataArrayInt::New(); neighbors->alloc(neighborsIdx->back(),1);
+  {
+    const int *neighIdx(neighborsIdx->begin());
+    int *neigh(neighbors->getPointer()),nodeId(0);
+    for(std::vector< std::set<int> >::const_iterator it=st0.begin();it!=st0.end();it++,neighIdx++,nodeId++)
+      {
+        std::set<int> s(*it); s.erase(nodeId);
+        std::copy(s.begin(),s.end(),neigh+*neighIdx);
+      }
+  }
+}
+
 /*!
  * Converts specified cells to either polygons (if \a this is a 2D mesh) or
  * polyhedrons (if \a this is a 3D mesh). The cells to convert are specified by an
@@ -954,7 +1027,7 @@ void MEDCouplingUMesh::convertToPolyTypes(const int *cellIdsToConvertBg, const i
   int nbOfCells(getNumberOfCells());
   if(dim==2)
     {
-      const int *connIndex=_nodal_connec_index->getConstPointer();
+      const int *connIndex=_nodal_connec_index->begin();
       int *conn=_nodal_connec->getPointer();
       for(const int *iter=cellIdsToConvertBg;iter!=cellIdsToConvertEnd;iter++)
         {
@@ -1241,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
@@ -1631,7 +1706,7 @@ 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
@@ -1828,7 +1903,7 @@ 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
@@ -1917,19 +1992,19 @@ void MEDCouplingUMesh::setPartOfMySelf(const int *cellIdsBg, const int *cellIdsE
       oss << ", whereas other mesh dimension is set equal to " << otherOnSameCoordsThanThis.getMeshDimension() << " !";
       throw INTERP_KERNEL::Exception(oss.str());
     }
-  int nbOfCellsToModify=(int)std::distance(cellIdsBg,cellIdsEnd);
+  std::size_t nbOfCellsToModify(std::distance(cellIdsBg,cellIdsEnd));
   if(nbOfCellsToModify!=otherOnSameCoordsThanThis.getNumberOfCells())
     {
       std::ostringstream oss; oss << "MEDCouplingUMesh::setPartOfMySelf : cells ids length (" <<  nbOfCellsToModify << ") do not match the number of cells of other mesh (" << otherOnSameCoordsThanThis.getNumberOfCells() << ") !";
       throw INTERP_KERNEL::Exception(oss.str());
     }
-  int nbOfCells=getNumberOfCells();
-  bool easyAssign=true;
-  const int *connI=_nodal_connec_index->getConstPointer();
-  const int *connIOther=otherOnSameCoordsThanThis._nodal_connec_index->getConstPointer();
+  std::size_t nbOfCells(getNumberOfCells());
+  bool easyAssign(true);
+  const int *connI(_nodal_connec_index->begin());
+  const int *connIOther=otherOnSameCoordsThanThis._nodal_connec_index->begin();
   for(const int *it=cellIdsBg;it!=cellIdsEnd && easyAssign;it++,connIOther++)
     {
-      if(*it>=0 && *it<nbOfCells)
+      if(*it>=0 && *it<(int)nbOfCells)
         {
           easyAssign=(connIOther[1]-connIOther[0])==(connI[*it+1]-connI[*it]);
         }
@@ -1967,7 +2042,7 @@ void MEDCouplingUMesh::setPartOfMySelfSlice(int start, int end, int step, const
       throw INTERP_KERNEL::Exception(oss.str());
     }
   int nbOfCellsToModify=DataArray::GetNumberOfItemGivenBESRelative(start,end,step,"MEDCouplingUMesh::setPartOfMySelfSlice : ");
-  if(nbOfCellsToModify!=otherOnSameCoordsThanThis.getNumberOfCells())
+  if(nbOfCellsToModify!=(int)otherOnSameCoordsThanThis.getNumberOfCells())
     {
       std::ostringstream oss; oss << "MEDCouplingUMesh::setPartOfMySelfSlice : cells ids length (" <<  nbOfCellsToModify << ") do not match the number of cells of other mesh (" << otherOnSameCoordsThanThis.getNumberOfCells() << ") !";
       throw INTERP_KERNEL::Exception(oss.str());
@@ -2210,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
@@ -2296,7 +2371,7 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On
   DAInt neighIInit00(tmp11);
   // Neighbor information of the mesh WITH the crack (some neighbors are removed):
   DataArrayInt *idsTmp=0;
-  bool b=m01->areCellsIncludedIn(&otherDimM1OnSameCoords,2,idsTmp);
+  m01->areCellsIncludedIn(&otherDimM1OnSameCoords,2,idsTmp);
   DAInt ids(idsTmp);
   // In the neighbor information remove the connection between high dimension cells and its low level constituents which are part
   // of the frontier given in parameter (i.e. the cells of low dimension from the group delimiting the crack):
@@ -2323,7 +2398,7 @@ void MEDCouplingUMesh::findNodesToDuplicate(const MEDCouplingUMesh& otherDimM1On
       // Connex zone without the crack (to compute the next seed really)
       int dnu;
       DAInt connexCheck = MEDCouplingUMesh::ComputeSpreadZoneGraduallyFromSeed(&seed, &seed+1, neighInit00,neighIInit00, -1, dnu);
-      int cnt = 0;
+      std::size_t cnt(0);
       for (int * ptr = connexCheck->getPointer(); cnt < connexCheck->getNumberOfTuples(); ptr++, cnt++)
         hitCells->setIJ(*ptr,0,1);
       // Connex zone WITH the crack (to identify cells lying on either part of the crack)
@@ -2355,7 +2430,7 @@ 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 
  * 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...
  * 
@@ -2496,7 +2571,7 @@ void MEDCouplingUMesh::shiftNodeNumbersInConn(int delta)
  * 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 
  * 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...
  * 
@@ -2708,11 +2783,10 @@ DataArrayInt *MEDCouplingUMesh::getCellsInBoundingBox(const INTERP_KERNEL::Direc
  *  \return INTERP_KERNEL::NormalizedCellType - enumeration item describing the cell type.
  *  \throw If \a cellId is invalid. Valid range is [0, \a this->getNumberOfCells() ).
  */
-INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::getTypeOfCell(int cellId) const
+INTERP_KERNEL::NormalizedCellType MEDCouplingUMesh::getTypeOfCell(std::size_t cellId) const
 {
-  const int *ptI=_nodal_connec_index->getConstPointer();
-  const int *pt=_nodal_connec->getConstPointer();
-  if(cellId>=0 && cellId<(int)_nodal_connec_index->getNbOfElems()-1)
+  const int *ptI(_nodal_connec_index->begin()),*pt(_nodal_connec->begin());
+  if(cellId<_nodal_connec_index->getNbOfElems()-1)
     return (INTERP_KERNEL::NormalizedCellType) pt[ptI[cellId]];
   else
     {
@@ -2754,13 +2828,11 @@ DataArrayInt *MEDCouplingUMesh::giveCellsWithType(INTERP_KERNEL::NormalizedCellT
 /*!
  * Returns nb of cells having the geometric type \a type. No throw if no cells in \a this has the geometric type \a type.
  */
-int MEDCouplingUMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
+std::size_t MEDCouplingUMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType type) const
 {
-  const int *ptI=_nodal_connec_index->getConstPointer();
-  const int *pt=_nodal_connec->getConstPointer();
-  int nbOfCells=getNumberOfCells();
-  int ret=0;
-  for(int i=0;i<nbOfCells;i++)
+  const int *ptI(_nodal_connec_index->begin()),*pt(_nodal_connec->begin());
+  std::size_t nbOfCells(getNumberOfCells()),ret(0);
+  for(std::size_t i=0;i<nbOfCells;i++)
     if((INTERP_KERNEL::NormalizedCellType) pt[ptI[i]]==type)
       ret++;
   return ret;
@@ -2775,10 +2847,9 @@ int MEDCouplingUMesh::getNumberOfCellsWithType(INTERP_KERNEL::NormalizedCellType
  *         cleared before the appending.
  *  \throw If \a cellId is invalid. Valid range is [0, \a this->getNumberOfCells() ).
  */
-void MEDCouplingUMesh::getNodeIdsOfCell(int cellId, std::vector<int>& conn) const
+void MEDCouplingUMesh::getNodeIdsOfCell(std::size_t cellId, std::vector<int>& conn) const
 {
-  const int *ptI=_nodal_connec_index->getConstPointer();
-  const int *pt=_nodal_connec->getConstPointer();
+  const int *ptI(_nodal_connec_index->begin()),*pt(_nodal_connec->begin());
   for(const int *w=pt+ptI[cellId]+1;w!=pt+ptI[cellId+1];w++)
     if(*w>=0)
       conn.push_back(*w);
@@ -2871,7 +2942,7 @@ 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).
@@ -2933,7 +3004,7 @@ 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
@@ -2977,14 +3048,14 @@ void MEDCouplingUMesh::setConnectivity(DataArrayInt *conn, DataArrayInt *connInd
  * Copy constructor. If 'deepCopy' is false \a this is a shallow copy of other.
  * If 'deeCpy' is true all arrays (coordinates and connectivities) are deeply copied.
  */
-MEDCouplingUMesh::MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCopy):MEDCouplingPointSet(other,deepCopy),_mesh_dim(other._mesh_dim),
+MEDCouplingUMesh::MEDCouplingUMesh(const MEDCouplingUMesh& other, bool deepCpy):MEDCouplingPointSet(other,deepCpy),_mesh_dim(other._mesh_dim),
     _nodal_connec(0),_nodal_connec_index(0),
     _types(other._types)
 {
   if(other._nodal_connec)
-    _nodal_connec=other._nodal_connec->performCopyOrIncrRef(deepCopy);
+    _nodal_connec=other._nodal_connec->performCopyOrIncrRef(deepCpy);
   if(other._nodal_connec_index)
-    _nodal_connec_index=other._nodal_connec_index->performCopyOrIncrRef(deepCopy);
+    _nodal_connec_index=other._nodal_connec_index->performCopyOrIncrRef(deepCpy);
 }
 
 MEDCouplingUMesh::~MEDCouplingUMesh()
@@ -3010,7 +3081,7 @@ void MEDCouplingUMesh::computeTypes()
  *  \return int - the number of cells in \a this mesh.
  *  \throw If the nodal connectivity of cells is not defined.
  */
-int MEDCouplingUMesh::getNumberOfCells() const
+std::size_t MEDCouplingUMesh::getNumberOfCells() const
 { 
   if(_nodal_connec_index)
     return _nodal_connec_index->getNumberOfTuples()-1;
@@ -3628,7 +3699,12 @@ MCAuto<MEDCouplingUMesh> MEDCouplingUMesh::clipSingle3DCellByPlane(const double
   std::vector<int> cut3DCurve(mDesc1->getNumberOfCells(),-2);
   for(const int *it=cellIds1D->begin();it!=cellIds1D->end();it++)
     cut3DCurve[*it]=-1;
-  mDesc1->split3DCurveWithPlane(origin,vec,eps,cut3DCurve);
+  bool sameNbNodes;
+  {
+    int oldNbNodes(mDesc1->getNumberOfNodes());
+    mDesc1->split3DCurveWithPlane(origin,vec,eps,cut3DCurve);
+    sameNbNodes=(mDesc1->getNumberOfNodes()==oldNbNodes);
+  }
   std::vector< std::pair<int,int> > cut3DSurf(mDesc2->getNumberOfCells());
   AssemblyForSplitFrom3DCurve(cut3DCurve,nodes,mDesc2->getNodalConnectivity()->begin(),mDesc2->getNodalConnectivityIndex()->begin(),
                               mDesc1->getNodalConnectivity()->begin(),mDesc1->getNodalConnectivityIndex()->begin(),
@@ -3644,7 +3720,7 @@ MCAuto<MEDCouplingUMesh> MEDCouplingUMesh::clipSingle3DCellByPlane(const double
   std::vector<std::vector<int> > res;
   buildSubCellsFromCut(cut3DSurf,desc2->begin(),descIndx2->begin(),mDesc1->getCoords()->begin(),eps,res);
   std::size_t sz(res.size());
-  if(res.size()==mDesc1->getNumberOfCells())
+  if(res.size()==mDesc1->getNumberOfCells() && sameNbNodes)
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::clipSingle3DCellByPlane : cell is not clipped !");
   for(std::size_t i=0;i<sz;i++)
     {
@@ -3713,7 +3789,6 @@ MCAuto<MEDCouplingUMesh> MEDCouplingUMesh::clipSingle3DCellByPlane(const double
   conn2I->pushBackSilent(conn2->getNumberOfTuples());
   ret2->setConnectivity(conn2,conn2I,true);
   ret2->checkConsistencyLight();
-  ret2->writeVTK("ret2.vtu");
   ret2->orientCorrectlyPolyhedrons();
   return ret2;
 }
@@ -3772,7 +3847,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
 {
@@ -3811,7 +3886,7 @@ void MEDCouplingUMesh::project1D(const double *pt, const double *v, double eps,
   MCAuto<MEDCouplingFieldDouble> f=buildDirectionVectorField();
   const double *fPtr=f->getArray()->getConstPointer();
   double tmp[3];
-  for(int i=0;i<getNumberOfCells();i++)
+  for(std::size_t i=0;i<getNumberOfCells();i++)
     {
       const double *tmp1=fPtr+3*i;
       tmp[0]=tmp1[1]*v[2]-tmp1[2]*v[1];
@@ -3901,7 +3976,7 @@ DataArrayDouble *MEDCouplingUMesh::distanceToPoints(const DataArrayDouble *pts,
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoints works only for spaceDim=meshDim+1 !");
   if(meshDim!=2 && meshDim!=1)
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::distanceToPoints : only mesh dimension 2 and 1 are implemented !");
-  if(pts->getNumberOfComponents()!=spaceDim)
+  if((int)pts->getNumberOfComponents()!=spaceDim)
     {
       std::ostringstream oss; oss << "MEDCouplingUMesh::distanceToPoints : input pts DataArrayDouble has " << pts->getNumberOfComponents() << " components whereas it should be equal to " << spaceDim << " (mesh spaceDimension) !";
       throw INTERP_KERNEL::Exception(oss.str());
@@ -4784,6 +4859,27 @@ void MEDCouplingUMesh::orientCorrectlyPolyhedrons()
   updateTime();
 }
 
+/*!
+ * 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.
+ */
+void MEDCouplingUMesh::invertOrientationOfAllCells()
+{
+  checkConnectivityFullyDefined();
+  std::set<INTERP_KERNEL::NormalizedCellType> gts(getAllGeoTypes());
+  int *conn(_nodal_connec->getPointer());
+  const int *conni(_nodal_connec_index->begin());
+  for(std::set<INTERP_KERNEL::NormalizedCellType>::const_iterator gt=gts.begin();gt!=gts.end();gt++)
+    {
+      INTERP_KERNEL::AutoCppPtr<INTERP_KERNEL::OrientationInverter> oi(INTERP_KERNEL::OrientationInverter::BuildInstanceFrom(*gt));
+      MCAuto<DataArrayInt> cwt(giveCellsWithType(*gt));
+      for(const int *it=cwt->begin();it!=cwt->end();it++)
+        oi->operate(conn+conni[*it]+1,conn+conni[*it+1]);
+    }
+  updateTime();
+}
+
 /*!
  * Finds and fixes incorrectly oriented linear extruded volumes (INTERP_KERNEL::NORM_HEXA8,
  * INTERP_KERNEL::NORM_PENTA6, INTERP_KERNEL::NORM_HEXGP12 etc) to respect the MED convention
@@ -5329,6 +5425,8 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTreeFast() const
 DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arcDetEps) const
 {
   checkFullyDefined();
+  INTERP_KERNEL::QuadraticPlanarArcDetectionPrecision arcPrec(arcDetEps);
+
   int spaceDim(getSpaceDimension()),mDim(getMeshDimension()),nbOfCells(getNumberOfCells());
   if(spaceDim!=2 || mDim!=2)
     throw INTERP_KERNEL::Exception("MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic : This method should be applied on mesh with mesh dimension equal to 2 and space dimension also equal to 2!");
@@ -5340,7 +5438,6 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree2DQuadratic(double arc
     {
       const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[*connI]));
       int sz(connI[1]-connI[0]-1);
-      INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=arcDetEps;
       std::vector<INTERP_KERNEL::Node *> nodes(sz);
       INTERP_KERNEL::QuadraticPolygon *pol(0);
       for(int j=0;j<sz;j++)
@@ -5377,6 +5474,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);
   MCAuto<DataArrayDouble> ret(DataArrayDouble::New()); ret->alloc(nbOfCells,2*spaceDim);
   double *bbox(ret->getPointer());
   const double *coords(_coords->begin());
@@ -5385,7 +5483,6 @@ DataArrayDouble *MEDCouplingUMesh::getBoundingBoxForBBTree1DQuadratic(double arc
     {
       const INTERP_KERNEL::CellModel& cm(INTERP_KERNEL::CellModel::GetCellModel((INTERP_KERNEL::NormalizedCellType)conn[*connI]));
       int sz(connI[1]-connI[0]-1);
-      INTERP_KERNEL::QUADRATIC_PLANAR::_arc_detection_precision=arcDetEps;
       std::vector<INTERP_KERNEL::Node *> nodes(sz);
       INTERP_KERNEL::Edge *edge(0);
       for(int j=0;j<sz;j++)
@@ -5436,7 +5533,7 @@ namespace MEDCouplingImpl
  * \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.
+ * This parameter is kept only for compatibility with other method listed above.
  */
 std::vector<int> MEDCouplingUMesh::getDistributionOfTypes() const
 {
@@ -5481,7 +5578,7 @@ std::vector<int> MEDCouplingUMesh::getDistributionOfTypes() const
  * 
  * 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.
+ * 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
 {
@@ -5584,7 +5681,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 !");
@@ -5626,7 +5723,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);
@@ -5792,7 +5889,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
 {
@@ -5831,7 +5928,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.
  */
@@ -5841,7 +5938,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.
@@ -5860,7 +5957,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()
 {
@@ -6033,7 +6130,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.
  *
@@ -6387,8 +6484,15 @@ DataArrayDouble *MEDCouplingUMesh::computePlaneEquationOf3DFaces() const
   for(int i=0;i<nbOfCells;i++,nodalI++,retPtr+=4)
     {
       double matrix[16]={0,0,0,1,0,0,0,1,0,0,0,1,1,1,1,0},matrix2[16];
-      if(nodalI[1]-nodalI[0]>=3)
+      if(nodalI[1]-nodalI[0]>=4)
         {
+          double aa[3]={coor[nodal[nodalI[0]+1+1]*3+0]-coor[nodal[nodalI[0]+1+0]*3+0],
+                        coor[nodal[nodalI[0]+1+1]*3+1]-coor[nodal[nodalI[0]+1+0]*3+1],
+                        coor[nodal[nodalI[0]+1+1]*3+2]-coor[nodal[nodalI[0]+1+0]*3+2]}
+          ,bb[3]={coor[nodal[nodalI[0]+1+2]*3+0]-coor[nodal[nodalI[0]+1+0]*3+0],
+                        coor[nodal[nodalI[0]+1+2]*3+1]-coor[nodal[nodalI[0]+1+0]*3+1],
+                        coor[nodal[nodalI[0]+1+2]*3+2]-coor[nodal[nodalI[0]+1+0]*3+2]};
+          double cc[3]={aa[1]*bb[2]-aa[2]*bb[1],aa[2]*bb[0]-aa[0]*bb[2],aa[0]*bb[1]-aa[1]*bb[0]};
           for(int j=0;j<3;j++)
             {
               int nodeId(nodal[nodalI[0]+1+j]);
@@ -6400,14 +6504,34 @@ DataArrayDouble *MEDCouplingUMesh::computePlaneEquationOf3DFaces() const
                   throw INTERP_KERNEL::Exception(oss.str());
                 }
             }
+          if(sqrt(cc[0]*cc[0]+cc[1]*cc[1]+cc[2]*cc[2])>1e-7)
+            {
+              INTERP_KERNEL::inverseMatrix(matrix,4,matrix2);
+              retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15];
+            }
+          else
+            {
+              if(nodalI[1]-nodalI[0]==4)
+                {
+                  std::ostringstream oss; oss << "MEDCouplingUMesh::computePlaneEquationOf3DFaces : cell" << i << " : Presence of The 3 colinear points !";
+                  throw INTERP_KERNEL::Exception(oss.str());
+                }
+              //
+              double dd[3]={0.,0.,0.};
+              for(int offset=nodalI[0]+1;offset<nodalI[1];offset++)
+                std::transform(coor+3*nodal[offset],coor+3*(nodal[offset]+1),dd,dd,std::plus<double>());
+              int nbOfNodesInCell(nodalI[1]-nodalI[0]-1);
+              std::transform(dd,dd+3,dd,std::bind2nd(std::multiplies<double>(),1./(double)nbOfNodesInCell));
+              std::copy(dd,dd+3,matrix+4*2);
+              INTERP_KERNEL::inverseMatrix(matrix,4,matrix2);
+              retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15];
+            }
         }
       else
         {
           std::ostringstream oss; oss << "MEDCouplingUMesh::computePlaneEquationOf3DFaces : invalid 2D cell #" << i << " ! Must be constitued by more than 3 nodes !";
           throw INTERP_KERNEL::Exception(oss.str());
         }
-      INTERP_KERNEL::inverseMatrix(matrix,4,matrix2);
-      retPtr[0]=matrix2[3]; retPtr[1]=matrix2[7]; retPtr[2]=matrix2[11]; retPtr[3]=matrix2[15];
     }
   return ret.retn();
 }
@@ -6669,7 +6793,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.
@@ -6723,7 +6847,7 @@ 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. 
  *  \param [in,out] meshes - a vector of meshes to update.
@@ -6774,7 +6898,7 @@ void MEDCouplingUMesh::MergeNodesOnUMeshesSharingSameCoords(const std::vector<ME
 
 
 /*!
- * 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)
 {
@@ -6818,7 +6942,7 @@ bool MEDCouplingUMesh::IsPolygonWellOriented(bool isQuadratic, const double *vec
 }
 
 /*!
- * 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)
 {
@@ -6904,48 +7028,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
@@ -6953,13 +7091,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)
@@ -7175,7 +7312,7 @@ DataArrayInt *MEDCouplingUMesh::buildUnionOf3DMesh() const
  *  means 6 arcs (0,1), (0,2), (0,3), (1,2), (1,3), (2,3)
  *  Arcs are not doubled but reflexive (1,1) arcs are present for each cell
  */
-MEDCouplingSkyLineArray *MEDCouplingUMesh::generateGraph() const
+MEDCouplingSkyLineArrayMEDCouplingUMesh::generateGraph() const
 {
   checkConnectivityFullyDefined();
 
@@ -7237,7 +7374,7 @@ MEDCouplingSkyLineArray *MEDCouplingUMesh::generateGraph() const
     cell2cell_index[icell+1]=cell2cell_index[icell]+cell2cell_index[icell+1];
 
   //filling up index and value to create skylinearray structure
-  MEDCouplingSkyLineArray* array=new MEDCouplingSkyLineArray(cell2cell_index,cell2cell);
+  MEDCouplingSkyLineArray * array(MEDCouplingSkyLineArray::New(cell2cell_index,cell2cell));
   return array;
 }
 
@@ -7611,7 +7748,7 @@ bool MEDCouplingUMesh::RemoveIdsFromIndexedArrays(const int *idsToRemoveBg, cons
       previousArrI=*arrIPtr;
       *arrIPtr=(int)arrOut.size();
     }
-  if(arr->getNumberOfTuples()==(int)arrOut.size())
+  if(arr->getNumberOfTuples()==arrOut.size())
     return false;
   arr->alloc((int)arrOut.size(),1);
   std::copy(arrOut.begin(),arrOut.end(),arr->getPointer());
@@ -7769,7 +7906,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)
@@ -7835,7 +7972,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)
@@ -7932,7 +8069,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)
@@ -7997,7 +8134,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)
@@ -8167,7 +8304,7 @@ DataArrayInt *MEDCouplingUMesh::ComputeRangesFromTypeDistribution(const std::vec
  *         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).