Salome HOME
Copyright update 2022
[tools/medcoupling.git] / src / INTERP_KERNEL / Geometric2D / InterpKernelGeo2DQuadraticPolygon.cxx
index a1df7ca5ce84c0c769c30c0d021e9af99b39c813..f6c1ee3164a5cbcf525d754ca07fad61c474a9b4 100644 (file)
@@ -1,4 +1,4 @@
-// Copyright (C) 2007-2014  CEA/DEN, EDF R&D
+// Copyright (C) 2007-2022  CEA/DEN, EDF R&D
 //
 // This library is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -47,7 +47,7 @@ QuadraticPolygon::QuadraticPolygon(const char *file)
   std::ifstream stream(file);
   stream.exceptions(std::ios_base::eofbit);
   try
-    {
+  {
       do
         stream.getline(currentLine,MAX_SIZE_OF_LINE_XFIG_FILE);
       while(strcmp(currentLine,"1200 2")!=0);
@@ -59,10 +59,25 @@ QuadraticPolygon::QuadraticPolygon(const char *file)
           pushBack(newEdge);
         }
       while(1);
-    }
-  catch(std::ifstream::failure&)
-    {
-    }
+  }
+  catch(const std::ifstream::failure&)
+  {
+  }
+  catch(const std::exception & ex)
+  {
+      // Some code before this catch throws the C++98 version of the exception (mangled
+      // name is " NSt8ios_base7failureE"), but FED24 compilation of the current version of the code
+      // tries to catch the C++11 version of it (mangled name "NSt8ios_base7failureB5cxx11E").
+      // So we have this nasty hack to catch both versions ...
+
+      // TODO: the below should be replaced by a better handling avoiding exception throwing.
+      if (std::string(ex.what()) == "basic_ios::clear")
+        {
+          //std::cout << "std::ios_base::failure C++11\n";
+        }
+      else
+        throw ex;
+  }
   front()->changeStartNodeWith(back()->getEndNode());
 }
 
@@ -252,11 +267,11 @@ double QuadraticPolygon::intersectWithAbs(QuadraticPolygon& other)
 }
 
 /*!
- * This method splits 'this' with 'other' into smaller pieces localizable. 'mapThis' is a map that gives the correspondance
+ * This method splits 'this' with 'other' into smaller pieces localizable. 'mapThis' is a map that gives the correspondence
  * between nodes contained in 'this' and node ids in a global mesh.
- * In the same way, 'mapOther' gives the correspondance between nodes contained in 'other' and node ids in a
- * global mesh from wich 'other' is extracted.
- * This method has 1 out paramater : 'edgesThis', After the call of this method, it contains the nodal connectivity (including type)
+ * In the same way, 'mapOther' gives the correspondence between nodes contained in 'other' and node ids in a
+ * global mesh from which 'other' is extracted.
+ * This method has 1 out parameter : 'edgesThis', After the call of this method, it contains the nodal connectivity (including type)
  * of 'this' into globlal "this mesh".
  * This method has 2 in/out parameters : 'subDivOther' and 'addCoo'.'otherEdgeIds' is useful to put values in
  * 'edgesThis', 'subDivOther' and 'addCoo'.
@@ -269,74 +284,96 @@ double QuadraticPolygon::intersectWithAbs(QuadraticPolygon& other)
  * the cell id in global other mesh.
  */
 void QuadraticPolygon::splitAbs(QuadraticPolygon& other,
-        const std::map<INTERP_KERNEL::Node *,int>& mapThis, const std::map<INTERP_KERNEL::Node *,int>& mapOther,
-        int offset1, int offset2 ,
-        const std::vector<int>& otherEdgeIds,
-        std::vector<int>& edgesThis, int cellIdThis,
-        std::vector< std::vector<int> >& edgesInOtherColinearWithThis, std::vector< std::vector<int> >& subDivOther,
-        std::vector<double>& addCoo)
+                                const std::map<INTERP_KERNEL::Node *,mcIdType>& mapThis, const std::map<INTERP_KERNEL::Node *,mcIdType>& mapOther,
+                                mcIdType offset1, mcIdType offset2 ,
+                                const std::vector<mcIdType>& otherEdgeIds,
+                                std::vector<mcIdType>& edgesThis, mcIdType cellIdThis,
+                                std::vector< std::vector<mcIdType> >& edgesInOtherColinearWithThis, std::vector< std::vector<mcIdType> >& subDivOther,
+                                std::vector<double>& addCoo, std::map<mcIdType,mcIdType>& mergedNodes)
 {
   double xBaryBB, yBaryBB;
   double fact=normalizeExt(&other, xBaryBB, yBaryBB);
+
   //
-  IteratorOnComposedEdge it1(this),it3(&other);
+  IteratorOnComposedEdge itThis(this),itOther(&other);  // other is (part of) the tool mesh
   MergePoints merge;
-  ComposedEdge *c1=new ComposedEdge;
-  ComposedEdge *c2=new ComposedEdge;
+  ComposedEdge *cThis=new ComposedEdge;
+  ComposedEdge *cOther=new ComposedEdge;
   int i=0;
-  std::map<INTERP_KERNEL::Node *,int> mapAddCoo;
-  for(it3.first();!it3.finished();it3.next(),i++)//iteration over 'other' _sub_edges
+  std::map<INTERP_KERNEL::Node *,mcIdType> mapAddCoo;
+  for(itOther.first();!itOther.finished();itOther.next(),i++)
     {
+      // For each edge of 'other', proceed with intersections: the edge might split into sub-edges, 'otherTmp' will hold the final split result.
+      // In the process of going through all edges of 'other', 'this' (which contains initially only one edge)
+      // is sub-divised into several edges : each of them has to be tested when intersecting the next candidate stored in 'other'.
       QuadraticPolygon otherTmp;
-      ElementaryEdge* curE3=it3.current();
-      otherTmp.pushBack(new ElementaryEdge(curE3->getPtr(),curE3->getDirection())); curE3->getPtr()->incrRef();
-      IteratorOnComposedEdge it2(&otherTmp);
-      for(it2.first();!it2.finished();it2.next())//iteration on subedges of 'other->_sub_edge'
+      ElementaryEdge* curOther=itOther.current();
+      otherTmp.pushBack(new ElementaryEdge(curOther->getPtr(),curOther->getDirection())); curOther->getPtr()->incrRef();
+      IteratorOnComposedEdge itOtherTmp(&otherTmp);
+      for(itOtherTmp.first();!itOtherTmp.finished();itOtherTmp.next())
         {
-          ElementaryEdge* curE2=it2.current();
-          if(!curE2->isThereStartPoint())
-            it1.first();
+          ElementaryEdge* curOtherTmp=itOtherTmp.current();
+          if(!curOtherTmp->isThereStartPoint())
+            itThis.first();   // reset iterator on 'this'
           else
-            it1=curE2->getIterator();
-          for(;!it1.finished();)//iteration over 'this' _sub_edges
+            itThis=curOtherTmp->getIterator();
+          for(;!itThis.finished();)
             {
-              ElementaryEdge* curE1=it1.current();
+              ElementaryEdge* curThis=itThis.current();
               merge.clear();
-              if(curE1->getPtr()->intersectWith(curE2->getPtr(),merge,*c1,*c2))
+              //
+              std::map<INTERP_KERNEL::Node *,mcIdType>::const_iterator thisStart(mapThis.find(curThis->getStartNode())),thisEnd(mapThis.find(curThis->getEndNode())),
+                                                                  otherStart(mapOther.find(curOtherTmp->getStartNode())),otherEnd(mapOther.find(curOtherTmp->getEndNode()));
+              mcIdType thisStart2(thisStart==mapThis.end()?-1:(*thisStart).second), thisEnd2(thisEnd==mapThis.end()?-1:(*thisEnd).second),
+                  otherStart2(otherStart==mapOther.end()?-1:(*otherStart).second+offset1),otherEnd2(otherEnd==mapOther.end()?-1:(*otherEnd).second+offset1);
+              //
+              if(curThis->getPtr()->intersectWith(curOtherTmp->getPtr(),merge,*cThis,*cOther))
                 {
-                  if(!curE1->getDirection()) c1->reverse();
-                  if(!curE2->getDirection()) c2->reverse();
-                  UpdateNeighbours(merge,it1,it2,c1,c2);
-                  //Substitution of simple edge by sub-edges.
-                  delete curE1; // <-- destroying simple edge coming from pol1
-                  delete curE2; // <-- destroying simple edge coming from pol2
-                  it1.insertElemEdges(c1,true);// <-- 2nd param is true to go next.
-                  it2.insertElemEdges(c2,false);// <-- 2nd param is false to avoid to go next.
-                  curE2=it2.current();
+                  if(!curThis->getDirection()) cThis->reverse();
+                  if(!curOtherTmp->getDirection()) cOther->reverse();
+                  // Substitution of a single simple edge by two sub-edges resulting from the intersection
+                  //    First modify the edges currently pointed by itThis and itOtherTmp so that the newly created node
+                  //    becomes the end of the previous sub-edge and the beginning of the next one.
+                  UpdateNeighbours(merge,itThis,itOtherTmp,cThis,cOther);
+                  delete curThis;           // <-- destroying simple edge coming from pol1
+                  delete curOtherTmp;       // <-- destroying simple edge coming from pol2
+                  //     Then insert second part of the intersection.
+                  itThis.insertElemEdges(cThis,true);       // <-- 2nd param is true to go next.
+                  itOtherTmp.insertElemEdges(cOther,false); // <-- 2nd param is false to avoid to go next.
+                  curOtherTmp=itOtherTmp.current();
                   //
-                  it1.assignMySelfToAllElems(c2);//To avoid that others
-                  SoftDelete(c1);
-                  SoftDelete(c2);
-                  c1=new ComposedEdge;
-                  c2=new ComposedEdge;
+                  itThis.assignMySelfToAllElems(cOther);
+                  SoftDelete(cThis);
+                  SoftDelete(cOther);
+                  cThis=new ComposedEdge;
+                  cOther=new ComposedEdge;
                 }
               else
                 {
-                  UpdateNeighbours(merge,it1,it2,curE1,curE2);
-                  it1.next();
+                  UpdateNeighbours(merge,itThis,itOtherTmp,curThis,curOtherTmp);
+                  itThis.next();
                 }
+              merge.updateMergedNodes(thisStart2,thisEnd2,otherStart2,otherEnd2,mergedNodes);
             }
         }
+      // If one sub-edge of otherTmp is "ON" an edge of this, then we have colinearity (all edges in otherTmp are //)
       if(otherTmp.presenceOfOn())
         edgesInOtherColinearWithThis[otherEdgeIds[i]].push_back(cellIdThis);
-      if(otherTmp._sub_edges.size()>1)
+      // Converting back to integer connectivity:
+      if(otherTmp._sub_edges.size()>1)   // only if a new point has been added (i.e. an actual intersection was done)
         {
-          for(std::list<ElementaryEdge *>::const_iterator it=otherTmp._sub_edges.begin();it!=otherTmp._sub_edges.end();it++)
-            (*it)->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,/**/fact,xBaryBB,yBaryBB,/**/subDivOther[otherEdgeIds[i]],addCoo,mapAddCoo);
+          std::size_t jj = 0, sz(otherTmp._sub_edges.size());
+          for(std::list<ElementaryEdge *>::const_iterator it=otherTmp._sub_edges.begin();it!=otherTmp._sub_edges.end();it++, jj++)
+            {
+              short skipStartOrEnd = jj == 0 ? -1 : (jj == sz-1 ? 1 : 0);  // -1 means START, 1 means END, 0 other
+              (*it)->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,
+                                      fact,xBaryBB,yBaryBB, skipStartOrEnd,
+                                      /*out*/ subDivOther[otherEdgeIds[i]],addCoo,mapAddCoo);
+            }
         }
     }
-  Delete(c1);
-  Delete(c2);
+  Delete(cThis);
+  Delete(cOther);
   //
   for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++)
     (*it)->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,/**/fact,xBaryBB,yBaryBB,/**/edgesThis,addCoo,mapAddCoo);
@@ -349,8 +386,8 @@ void QuadraticPolygon::splitAbs(QuadraticPolygon& other,
  * orientation of edge (see buildDescendingConnectivity2() method).
  * See appendEdgeFromCrudeDataArray() for params description.
  */
-void QuadraticPolygon::buildFromCrudeDataArray(const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad, const int *nodalBg, const double *coords,
-                                               const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges)
+void QuadraticPolygon::buildFromCrudeDataArray(const std::map<mcIdType,INTERP_KERNEL::Node *>& mapp, bool isQuad, const mcIdType *nodalBg, const double *coords,
+                                               const mcIdType *descBg, const mcIdType *descEnd, const std::vector<std::vector<mcIdType> >& intersectEdges)
 {
   std::size_t nbOfSeg=std::distance(descBg,descEnd);
   for(std::size_t i=0;i<nbOfSeg;i++)
@@ -359,15 +396,15 @@ void QuadraticPolygon::buildFromCrudeDataArray(const std::map<int,INTERP_KERNEL:
     }
 }
 
-void QuadraticPolygon::appendEdgeFromCrudeDataArray(std::size_t edgePos, const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad,
-                                                    const int *nodalBg, const double *coords,
-                                                    const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges)
+void QuadraticPolygon::appendEdgeFromCrudeDataArray(std::size_t edgePos, const std::map<mcIdType,INTERP_KERNEL::Node *>& mapp, bool isQuad,
+                                                    const mcIdType *nodalBg, const double *coords,
+                                                    const mcIdType *descBg, const mcIdType *descEnd, const std::vector<std::vector<mcIdType> >& intersectEdges)
 {
   if(!isQuad)
     {
       bool direct=descBg[edgePos]>0;
-      int edgeId=abs(descBg[edgePos])-1; // back to C indexing mode
-      const std::vector<int>& subEdge=intersectEdges[edgeId];
+      mcIdType edgeId=std::abs(descBg[edgePos])-1; // back to C indexing mode
+      const std::vector<mcIdType>& subEdge=intersectEdges[edgeId];
       std::size_t nbOfSubEdges=subEdge.size()/2;
       for(std::size_t j=0;j<nbOfSubEdges;j++)
         appendSubEdgeFromCrudeDataArray(0,j,direct,edgeId,subEdge,mapp);
@@ -389,8 +426,8 @@ void QuadraticPolygon::appendEdgeFromCrudeDataArray(std::size_t edgePos, const s
       delete e1; delete e2;
       //
       bool direct=descBg[edgePos]>0;
-      int edgeId=abs(descBg[edgePos])-1;
-      const std::vector<int>& subEdge=intersectEdges[edgeId];
+      mcIdType edgeId=std::abs(descBg[edgePos])-1;
+      const std::vector<mcIdType>& subEdge=intersectEdges[edgeId];
       std::size_t nbOfSubEdges=subEdge.size()/2;
       if(colinearity)
         {   
@@ -408,14 +445,14 @@ void QuadraticPolygon::appendEdgeFromCrudeDataArray(std::size_t edgePos, const s
     }
 }
 
-void QuadraticPolygon::appendSubEdgeFromCrudeDataArray(Edge *baseEdge, std::size_t j, bool direct, int edgeId, const std::vector<int>& subEdge, const std::map<int,INTERP_KERNEL::Node *>& mapp)
+void QuadraticPolygon::appendSubEdgeFromCrudeDataArray(Edge *baseEdge, std::size_t j, bool direct, mcIdType edgeId, const std::vector<mcIdType>& subEdge, const std::map<mcIdType,INTERP_KERNEL::Node *>& mapp)
 {
   std::size_t nbOfSubEdges=subEdge.size()/2;
   if(!baseEdge)
     {//it is not a quadratic subedge
       Node *start=(*mapp.find(direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1])).second;
       Node *end=(*mapp.find(direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2])).second;
-      ElementaryEdge *e=ElementaryEdge::BuildEdgeFromCrudeDataArray(true,start,end);
+      ElementaryEdge *e=ElementaryEdge::BuildEdgeFromStartEndDir(true,start,end);
       pushBack(e);
     }
   else
@@ -432,17 +469,17 @@ void QuadraticPolygon::appendSubEdgeFromCrudeDataArray(Edge *baseEdge, std::size
  * This method builds from descending conn of a quadratic polygon stored in crude mode (MEDCoupling). Descending conn is in FORTRAN relative mode in order to give the
  * orientation of edge.
  */
-void QuadraticPolygon::buildFromCrudeDataArray2(const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad, const int *nodalBg, const double *coords, const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges,
-                                                const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, const std::vector<std::vector<int> >& intersectEdges1,
-                                                const std::vector< std::vector<int> >& colinear1,
-                                                std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> >& alreadyExistingIn2)
+void QuadraticPolygon::buildFromCrudeDataArray2(const std::map<mcIdType,INTERP_KERNEL::Node *>& mapp, bool isQuad, const mcIdType *nodalBg, const double *coords, const mcIdType *descBg, const mcIdType *descEnd, const std::vector<std::vector<mcIdType> >& intersectEdges2,
+                                                const INTERP_KERNEL::QuadraticPolygon& pol1, const mcIdType *descBg1, const mcIdType *descEnd1, const std::vector<std::vector<mcIdType> >& intersectEdges1,
+                                                const std::vector< std::vector<mcIdType> >& colinear1,
+                                                std::map<mcIdType,std::vector<INTERP_KERNEL::ElementaryEdge *> >& alreadyExistingIn2)
 {
   std::size_t nbOfSeg=std::distance(descBg,descEnd);
   for(std::size_t i=0;i<nbOfSeg;i++)//loop over all edges of pol2
     {
       bool direct=descBg[i]>0;
-      int edgeId=abs(descBg[i])-1;//current edge id of pol2
-      std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> >::const_iterator it1=alreadyExistingIn2.find(descBg[i]),it2=alreadyExistingIn2.find(-descBg[i]);
+      mcIdType edgeId=std::abs(descBg[i])-1;//current edge id of pol2
+      std::map<mcIdType,std::vector<INTERP_KERNEL::ElementaryEdge *> >::const_iterator it1=alreadyExistingIn2.find(descBg[i]),it2=alreadyExistingIn2.find(-descBg[i]);
       if(it1!=alreadyExistingIn2.end() || it2!=alreadyExistingIn2.end())
         {
           bool sameDir=(it1!=alreadyExistingIn2.end());
@@ -466,28 +503,28 @@ void QuadraticPolygon::buildFromCrudeDataArray2(const std::map<int,INTERP_KERNEL
           continue;
         }
       bool directos=colinear1[edgeId].empty();
-      std::vector<std::pair<int,std::pair<bool,int> > > idIns1;
-      int offset1=0;
+      std::vector<std::pair<mcIdType,std::pair<bool,mcIdType> > > idIns1;
+      mcIdType offset1=0;
       if(!directos)
         {// if the current edge of pol2 has one or more colinear edges part into pol1
-          const std::vector<int>& c=colinear1[edgeId];
+          const std::vector<mcIdType>& c=colinear1[edgeId];
           std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1);
           for(std::size_t j=0;j<nbOfEdgesIn1;j++)
             {
-              int edgeId1=abs(descBg1[j])-1;
+              mcIdType edgeId1=std::abs(descBg1[j])-1;
               if(std::find(c.begin(),c.end(),edgeId1)!=c.end())
                 {
-                  idIns1.push_back(std::pair<int,std::pair<bool,int> >(edgeId1,std::pair<bool,int>(descBg1[j]>0,offset1)));// it exists an edge into pol1 given by tuple (idIn1,direct1) that is colinear at edge 'edgeId' in pol2
+                  idIns1.push_back(std::pair<mcIdType,std::pair<bool,mcIdType> >(edgeId1,std::pair<bool,mcIdType>(descBg1[j]>0,offset1)));// it exists an edge into pol1 given by tuple (idIn1,direct1) that is colinear at edge 'edgeId' in pol2
                   //std::pair<edgeId1); direct1=descBg1[j]>0;
                 }
-              offset1+=intersectEdges1[edgeId1].size()/2;//offset1 is used to find the INTERP_KERNEL::Edge * instance into pol1 that will be part of edge into pol2
+              offset1+=ToIdType(intersectEdges1[edgeId1].size()/2);//offset1 is used to find the INTERP_KERNEL::Edge * instance into pol1 that will be part of edge into pol2
             }
           directos=idIns1.empty();
         }
       if(directos)
         {//no subpart of edge 'edgeId' of pol2 is in pol1 so let's operate the same thing that QuadraticPolygon::buildFromCrudeDataArray method
           std::size_t oldSz=_sub_edges.size();
-          appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
+          appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges2);
           std::size_t newSz=_sub_edges.size();
           std::size_t zeSz=newSz-oldSz;
           alreadyExistingIn2[descBg[i]].resize(zeSz);
@@ -497,25 +534,25 @@ void QuadraticPolygon::buildFromCrudeDataArray2(const std::map<int,INTERP_KERNEL
         }
       else
         {//there is subpart of edge 'edgeId' of pol2 inside pol1
-          const std::vector<int>& subEdge=intersectEdges[edgeId];
+          const std::vector<mcIdType>& subEdge=intersectEdges2[edgeId];
           std::size_t nbOfSubEdges=subEdge.size()/2;
           for(std::size_t j=0;j<nbOfSubEdges;j++)
             {
-              int idBg=direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1];
-              int idEnd=direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2];
-              bool direction11,found=false;
-              bool direct1;//store if needed the direction in 1
-              int offset2;
-              std::size_t nbOfSubEdges1;
-              for(std::vector<std::pair<int,std::pair<bool,int> > >::const_iterator it=idIns1.begin();it!=idIns1.end() && !found;it++)
+              mcIdType idBg=direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1];
+              mcIdType idEnd=direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2];
+              bool direction11=false,found=false;
+              bool direct1=false;//store if needed the direction in 1
+              mcIdType offset2=0;
+              mcIdType nbOfSubEdges1=0;
+              for(std::vector<std::pair<mcIdType,std::pair<bool,mcIdType> > >::const_iterator it=idIns1.begin();it!=idIns1.end() && !found;it++)
                 {
-                  int idIn1=(*it).first;//store if needed the cell id in 1
+                  mcIdType idIn1=(*it).first;//store if needed the cell id in 1
                   direct1=(*it).second.first;
                   offset1=(*it).second.second;
-                  const std::vector<int>& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1];
-                  nbOfSubEdges1=subEdge1PossiblyAlreadyIn1.size()/2;
+                  const std::vector<mcIdType>& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1];
+                  nbOfSubEdges1=ToIdType(subEdge1PossiblyAlreadyIn1.size()/2);
                   offset2=0;
-                  for(std::size_t k=0;k<nbOfSubEdges1 && !found;k++)
+                  for(mcIdType k=0;k<nbOfSubEdges1 && !found;k++)
                     {//perform a loop on all subedges of pol1 that includes edge 'edgeId' of pol2. For the moment we iterate only on subedges of ['idIn1']... To improve
                       if(subEdge1PossiblyAlreadyIn1[2*k]==idBg && subEdge1PossiblyAlreadyIn1[2*k+1]==idEnd)
                         { direction11=true; found=true; }
@@ -530,13 +567,13 @@ void QuadraticPolygon::buildFromCrudeDataArray2(const std::map<int,INTERP_KERNEL
                   //appendEdgeFromCrudeDataArray(j,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
                   Node *start=(*mapp.find(idBg)).second;
                   Node *end=(*mapp.find(idEnd)).second;
-                  ElementaryEdge *e=ElementaryEdge::BuildEdgeFromCrudeDataArray(true,start,end);
+                  ElementaryEdge *e=ElementaryEdge::BuildEdgeFromStartEndDir(true,start,end);
                   pushBack(e);
                   alreadyExistingIn2[descBg[i]].push_back(e);
                 }
               else
                 {//the current subedge of edge 'edgeId' of pol2 is part of the colinear edge 'idIn1' of pol1 -> reuse Edge instance of pol1
-                  ElementaryEdge *e=pol1[offset1+(direct1?offset2:nbOfSubEdges1-offset2-1)];
+                  ElementaryEdge *e=pol1[FromIdType<int>(offset1+(direct1?offset2:nbOfSubEdges1-offset2-1))];
                   Edge *ee=e->getPtr();
                   ee->incrRef();
                   ElementaryEdge *e2=new ElementaryEdge(ee,!(direct1^direction11));
@@ -552,39 +589,39 @@ void QuadraticPolygon::buildFromCrudeDataArray2(const std::map<int,INTERP_KERNEL
  * Method expected to be called on pol2. Every params not suffixed by numbered are supposed to refer to pol2 (this).
  * Method to find edges that are ON.
  */
-void QuadraticPolygon::updateLocOfEdgeFromCrudeDataArray2(const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges,
-      const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1,
-      const std::vector<std::vector<int> >& intersectEdges1, const std::vector< std::vector<int> >& colinear1) const
+void QuadraticPolygon::updateLocOfEdgeFromCrudeDataArray2(const mcIdType *descBg, const mcIdType *descEnd, const std::vector<std::vector<mcIdType> >& intersectEdges,
+                                                          const INTERP_KERNEL::QuadraticPolygon& pol1, const mcIdType *descBg1, const mcIdType *descEnd1,
+                                                          const std::vector<std::vector<mcIdType> >& intersectEdges1, const std::vector< std::vector<mcIdType> >& colinear1) const
 {
   std::size_t nbOfSeg=std::distance(descBg,descEnd);
   for(std::size_t i=0;i<nbOfSeg;i++)//loop over all edges of pol2
     {
       bool direct=descBg[i]>0;
-      int edgeId=abs(descBg[i])-1;//current edge id of pol2
-      const std::vector<int>& c=colinear1[edgeId];
+      mcIdType edgeId=std::abs(descBg[i])-1;//current edge id of pol2
+      const std::vector<mcIdType>& c=colinear1[edgeId];
       if(c.empty())
         continue;
-      const std::vector<int>& subEdge=intersectEdges[edgeId];
+      const std::vector<mcIdType>& subEdge=intersectEdges[edgeId];
       std::size_t nbOfSubEdges=subEdge.size()/2;
       //
       std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1);
-      int offset1=0;
+      mcIdType offset1=0;
       for(std::size_t j=0;j<nbOfEdgesIn1;j++)
         {
-          int edgeId1=abs(descBg1[j])-1;
+          mcIdType edgeId1=std::abs(descBg1[j])-1;
           if(std::find(c.begin(),c.end(),edgeId1)!=c.end())
             {
               for(std::size_t k=0;k<nbOfSubEdges;k++)
                 {
-                  int idBg=direct?subEdge[2*k]:subEdge[2*nbOfSubEdges-2*k-1];
-                  int idEnd=direct?subEdge[2*k+1]:subEdge[2*nbOfSubEdges-2*k-2];
-                  int idIn1=edgeId1;
+                  mcIdType idBg=direct?subEdge[2*k]:subEdge[2*nbOfSubEdges-2*k-1];
+                  mcIdType idEnd=direct?subEdge[2*k+1]:subEdge[2*nbOfSubEdges-2*k-2];
+                  mcIdType idIn1=edgeId1;
                   bool direct1=descBg1[j]>0;
-                  const std::vector<int>& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1];
-                  std::size_t nbOfSubEdges1=subEdge1PossiblyAlreadyIn1.size()/2;
-                  int offset2=0;
+                  const std::vector<mcIdType>& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1];
+                  mcIdType nbOfSubEdges1=ToIdType(subEdge1PossiblyAlreadyIn1.size()/2);
+                  mcIdType offset2=0;
                   bool found=false;
-                  for(std::size_t kk=0;kk<nbOfSubEdges1 && !found;kk++)
+                  for(mcIdType kk=0;kk<nbOfSubEdges1 && !found;kk++)
                     {
                       found=(subEdge1PossiblyAlreadyIn1[2*kk]==idBg && subEdge1PossiblyAlreadyIn1[2*kk+1]==idEnd) || (subEdge1PossiblyAlreadyIn1[2*kk]==idEnd && subEdge1PossiblyAlreadyIn1[2*kk+1]==idBg);
                       if(!found)
@@ -592,17 +629,17 @@ void QuadraticPolygon::updateLocOfEdgeFromCrudeDataArray2(const int *descBg, con
                     }
                   if(found)
                     {
-                      ElementaryEdge *e=pol1[offset1+(direct1?offset2:nbOfSubEdges1-offset2-1)];
+                      ElementaryEdge *e=pol1[(int)(offset1+(direct1?offset2:nbOfSubEdges1-offset2-1))];
                       e->getPtr()->declareOn();
                     }
                 }
             }
-          offset1+=intersectEdges1[edgeId1].size()/2;//offset1 is used to find the INTERP_KERNEL::Edge * instance into pol1 that will be part of edge into pol2
+          offset1+=ToIdType(intersectEdges1[edgeId1].size()/2);//offset1 is used to find the INTERP_KERNEL::Edge * instance into pol1 that will be part of edge into pol2
         }
     }
 }
 
-void QuadraticPolygon::appendCrudeData(const std::map<INTERP_KERNEL::Node *,int>& mapp, double xBary, double yBary, double fact, int offset, std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI) const
+void QuadraticPolygon::appendCrudeData(const std::map<INTERP_KERNEL::Node *,mcIdType>& mapp, double xBary, double yBary, double fact, mcIdType offset, std::vector<double>& addCoordsQuadratic, std::vector<mcIdType>& conn, std::vector<mcIdType>& connI) const
 {
   int nbOfNodesInPg=0;
   bool presenceOfQuadratic=presenceOfQuadraticEdge();
@@ -611,14 +648,14 @@ void QuadraticPolygon::appendCrudeData(const std::map<INTERP_KERNEL::Node *,int>
     {
       Node *tmp=0;
       tmp=(*it)->getStartNode();
-      std::map<INTERP_KERNEL::Node *,int>::const_iterator it1=mapp.find(tmp);
+      std::map<INTERP_KERNEL::Node *,mcIdType>::const_iterator it1=mapp.find(tmp);
       conn.push_back((*it1).second);
       nbOfNodesInPg++;
     }
   if(presenceOfQuadratic)
     {
       int j=0;
-      int off=offset+((int)addCoordsQuadratic.size())/2;
+      mcIdType off=offset+ToIdType(addCoordsQuadratic.size())/2;
       for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++,j++,nbOfNodesInPg++)
         {
           INTERP_KERNEL::Node *node=(*it)->getPtr()->buildRepresentantOfMySelf();
@@ -633,18 +670,21 @@ void QuadraticPolygon::appendCrudeData(const std::map<INTERP_KERNEL::Node *,int>
 }
 
 /*!
- * This method make the hypothesis that 'this' and 'other' are splited at the minimum into edges that are fully IN, OUT or ON.
- * This method returns newly created polygons in 'conn' and 'connI' and the corresponding ids ('idThis','idOther') are stored respectively into 'nbThis' and 'nbOther'.
- * @param [in,out] edgesThis, parameter that keep informed the caller abount the edges in this not shared by the result of intersection of \a this with \a other
- * @param [in,out] edgesBoundaryOther, parameter that strores all edges in result of intersection that are not 
+ * This method make the hypothesis that \a this and \a other are split at the minimum into edges that are fully IN, OUT or ON.
+ * This method returns newly created polygons in \a conn and \a connI and the corresponding ids ( \a idThis, \a idOther) are stored respectively into \a nbThis and \a nbOther.
+ * @param [in,out] edgesThis parameter that keep informed the caller about the edges in this not shared by the result of intersection of \a this with \a other
+ * @param [in,out] edgesBoundaryOther parameter that stores all edges in result of intersection that are not
  */
-void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, std::set<INTERP_KERNEL::Edge *>& edgesThis, std::set<INTERP_KERNEL::Edge *>& edgesBoundaryOther, const std::map<INTERP_KERNEL::Node *,int>& mapp, int idThis, int idOther, int offset, std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI, std::vector<int>& nbThis, std::vector<int>& nbOther)
+void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, std::set<INTERP_KERNEL::Edge *>& edgesThis, std::set<INTERP_KERNEL::Edge *>& edgesBoundaryOther,
+                                          const std::map<INTERP_KERNEL::Node *,mcIdType>& mapp, mcIdType idThis, mcIdType idOther, mcIdType offset,
+                                          std::vector<double>& addCoordsQuadratic, std::vector<mcIdType>& conn, std::vector<mcIdType>& connI,
+                                          std::vector<mcIdType>& nbThis, std::vector<mcIdType>& nbOther)
 {
   double xBaryBB, yBaryBB;
   double fact=normalizeExt(&other, xBaryBB, yBaryBB);
-  //Locate 'this' relative to 'other'
+  //Locate \a this relative to \a other (edges of \a this, aka \a pol1 are marked as IN or OUT)
   other.performLocatingOperationSlow(*this);  // without any assumption
-  std::vector<QuadraticPolygon *> res=buildIntersectionPolygons(other,*this);
+  std::vector<QuadraticPolygon *> res=buildIntersectionPolygons(*this,other);
   for(std::vector<QuadraticPolygon *>::iterator it=res.begin();it!=res.end();it++)
     {
       (*it)->appendCrudeData(mapp,xBaryBB,yBaryBB,fact,offset,addCoordsQuadratic,conn,connI);
@@ -669,6 +709,28 @@ void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, std::set<INTE
   unApplyGlobalSimilarityExt(other,xBaryBB,yBaryBB,fact);
 }
 
+/*!
+ * Remove the two 'identical' edges from the list, and drecRef the objects.
+ */
+void QuadraticPolygon::cleanDegeneratedConsecutiveEdges()
+{
+  IteratorOnComposedEdge it(this);
+  ElementaryEdge * prevEdge = 0;
+  if  (recursiveSize() > 2)
+    for(it.first();!it.finished();it.next())
+      {
+        ElementaryEdge * cur = it.current();
+        if (prevEdge && prevEdge->hasSameExtremities(*cur))
+          {
+            it.eraseCurrent();
+            it.eraseCurrent();
+            prevEdge = it.current();
+          }
+        else
+            prevEdge = cur;
+      }
+}
+
 /*!
  * Warning This method is \b NOT const. 'this' and 'other' are modified after call of this method.
  * 'other' is a QuadraticPolygon of \b non closed edges.
@@ -688,7 +750,7 @@ double QuadraticPolygon::intersectWithAbs1D(QuadraticPolygon& other, bool& isCol
   for(std::list<ElementaryEdge *>::const_iterator it=cpyOfOther._sub_edges.begin();it!=cpyOfOther._sub_edges.end();it++)
     {
       switch((*it)->getLoc())
-        {
+      {
         case FULL_IN_1:
           {
             ret += fabs((*it)->getPtr()->getCurveLength());
@@ -703,7 +765,7 @@ double QuadraticPolygon::intersectWithAbs1D(QuadraticPolygon& other, bool& isCol
         default:
           {
           }
-        }
+      }
     }
   return ret * fact;
 }
@@ -730,7 +792,7 @@ double QuadraticPolygon::intersectWithAbs(QuadraticPolygon& other, double* baryc
     {
       barycenter[0]=barycenter[0]/ret*fact+xBaryBB;
       barycenter[1]=barycenter[1]/ret*fact+yBaryBB;
-      
+
     }
   return ret*fact*fact;
 }
@@ -868,6 +930,7 @@ void QuadraticPolygon::intersectForPoint(const QuadraticPolygon& other, std::vec
  * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
  * This is possible because loc attribute in Edge class is mutable.
  * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
+ * This method is currently not used by any high level functionality.
  */
 std::vector<QuadraticPolygon *> QuadraticPolygon::intersectMySelfWith(const QuadraticPolygon& other) const
 {
@@ -876,7 +939,7 @@ std::vector<QuadraticPolygon *> QuadraticPolygon::intersectMySelfWith(const Quad
   SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits);
   //At this point cpyOfThis and cpyOfOther have been splited at maximum edge so that in/out can been done.
   performLocatingOperation(cpyOfOther);
-  return other.buildIntersectionPolygons(cpyOfThis,cpyOfOther);
+  return other.buildIntersectionPolygons(cpyOfOther, cpyOfThis);
 }
 
 /*!
@@ -899,7 +962,7 @@ void QuadraticPolygon::SplitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticP
         it1=curE2->getIterator();
       for(;!it1.finished();)
         {
-          
+
           ElementaryEdge* curE1=it1.current();
           merge.clear(); nbOfSplits++;
           if(curE1->getPtr()->intersectWith(curE2->getPtr(),merge,*c1,*c2))
@@ -931,14 +994,14 @@ void QuadraticPolygon::SplitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticP
   Delete(c2);
 }
 
-void QuadraticPolygon::performLocatingOperation(QuadraticPolygon& pol2) const
+void QuadraticPolygon::performLocatingOperation(QuadraticPolygon& pol1) const
 {
-  IteratorOnComposedEdge it(&pol2);
+  IteratorOnComposedEdge it(&pol1);
   TypeOfEdgeLocInPolygon loc=FULL_ON_1;
   for(it.first();!it.finished();it.next())
     {
       ElementaryEdge *cur=it.current();
-      loc=cur->locateFullyMySelf(*this,loc);
+      loc=cur->locateFullyMySelf(*this,loc);//*this=pol2=other
     }
 }
 
@@ -953,34 +1016,36 @@ void QuadraticPolygon::performLocatingOperationSlow(QuadraticPolygon& pol2) cons
 }
 
 /*!
- * Given 2 polygons 'pol1' and 'pol2' (localized) the resulting polygons are returned.
+ * Given 2 polygons \a pol1 and \a pol2 (localized) the resulting polygons are returned.
  *
- * this : pol2 simplified.
- * @param pol1 pol1 split.
- * @param pol2 pol2 split.
+ * this : pol1 simplified.
+ * @param [in] pol1 pol1 split.
+ * @param [in] pol2 pol2 split.
  */
 std::vector<QuadraticPolygon *> QuadraticPolygon::buildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2) const
 {
   std::vector<QuadraticPolygon *> ret;
-  std::list<QuadraticPolygon *> pol2Zip=pol2.zipConsecutiveInSegments();
-  if(!pol2Zip.empty())
-    closePolygons(pol2Zip,pol1,ret);
+  // Extract from pol1, and pol1 only, all consecutive edges.
+  // pol1Zip contains concatenated pieces of pol1 which are part of the resulting intersecting cell being built.
+  std::list<QuadraticPolygon *> pol1Zip=pol1.zipConsecutiveInSegments();
+  if(!pol1Zip.empty())
+    ClosePolygons(pol1Zip,*this,pol2,ret);
   else
-    {//borders of pol2 do not cross pol1,and pol2 borders are outside of pol1. That is to say, either pol2 and pol1
-      //do not overlap or  pol1 is fully inside pol2. So in the first case no intersection, in the other case
-      //the intersection is pol1.
-      ElementaryEdge *e1FromPol1=pol1[0];
+    {//borders of pol1 do not cross pol2,and pol1 borders are outside of pol2. That is to say, either pol1 and pol2
+      //do not overlap or  pol2 is fully inside pol1. So in the first case no intersection, in the other case
+      //the intersection is pol2.
+      ElementaryEdge *e1FromPol2=pol2[0];
       TypeOfEdgeLocInPolygon loc=FULL_ON_1;
-      loc=e1FromPol1->locateFullyMySelf(*this,loc);
+      loc=e1FromPol2->locateFullyMySelf(*this,loc);
       if(loc==FULL_IN_1)
-        ret.push_back(new QuadraticPolygon(pol1));
+        ret.push_back(new QuadraticPolygon(pol2));
     }
   return ret;
 }
 
 /*!
  * Returns parts of potentially non closed-polygons. Each returned polygons are not mergeable.
- * this : pol2 split and locallized.
+ * this : pol1 split and localized.
  */
 std::list<QuadraticPolygon *> QuadraticPolygon::zipConsecutiveInSegments() const
 {
@@ -1015,129 +1080,147 @@ std::list<QuadraticPolygon *> QuadraticPolygon::zipConsecutiveInSegments() const
 }
 
 /*!
- * 'this' should be considered as pol2Simplified.
- * @param pol2zip is a list of set of edges (openned polygon) coming from split polygon 2.
- * @param pol1 is split pol1.
- * @param results the resulting \b CLOSED polygons.
+ * @param [in] pol1Zip is a list of set of edges (=an opened polygon) coming from split polygon 1.
+ * @param [in] pol1 should be considered as pol1Simplified.
+ * @param [in] pol2 is split pol2.
+ * @param [out] results the resulting \b CLOSED polygons.
  */
-void QuadraticPolygon::closePolygons(std::list<QuadraticPolygon *>& pol2Zip, const QuadraticPolygon& pol1,
-                                     std::vector<QuadraticPolygon *>& results) const
+void QuadraticPolygon::ClosePolygons(std::list<QuadraticPolygon *>& pol1Zip, const QuadraticPolygon& pol1, const QuadraticPolygon& pol2,
+                                     std::vector<QuadraticPolygon *>& results)
 {
-  bool directionKnownInPol1=false;
-  bool directionInPol1;
-  for(std::list<QuadraticPolygon *>::iterator iter=pol2Zip.begin();iter!=pol2Zip.end();)
+  bool directionKnownInPol2=false;
+  bool directionInPol2=false;
+  bool needCleaning = false;
+  for(std::list<QuadraticPolygon *>::iterator iter=pol1Zip.begin();iter!=pol1Zip.end();)
     {
+      // Build incrementally the full closed cells from the consecutive line parts already built in pol1Zip.
+      // At the end of the process the item currently iterated has been totally completed (start_node=end_node)
+      // This process can produce several cells.
       if((*iter)->completed())
         {
+          if (needCleaning)
+            (*iter)->cleanDegeneratedConsecutiveEdges();
           results.push_back(*iter);
-          directionKnownInPol1=false;
-          iter=pol2Zip.erase(iter);
+          directionKnownInPol2=false;
+          needCleaning=false;
+          iter=pol1Zip.erase(iter);
           continue;
         }
-      if(!directionKnownInPol1)
+      if(!directionKnownInPol2)
         {
-          if(!(*iter)->amIAChanceToBeCompletedBy(pol1,*this,directionInPol1))
-            { delete *iter; iter=pol2Zip.erase(iter); continue; }
+          if(!(*iter)->haveIAChanceToBeCompletedBy(pol1,pol2,directionInPol2, needCleaning))
+            { delete *iter; iter=pol1Zip.erase(iter); continue; }
           else
-            directionKnownInPol1=true;
+            directionKnownInPol2=true;
         }
       std::list<QuadraticPolygon *>::iterator iter2=iter; iter2++;
-      std::list<QuadraticPolygon *>::iterator iter3=(*iter)->fillAsMuchAsPossibleWith(pol1,iter2,pol2Zip.end(),directionInPol1);
-      if(iter3!=pol2Zip.end())
+      // Fill as much as possible the current iterate (=a part of pol1) with consecutive pieces from pol2:
+      std::list<QuadraticPolygon *>::iterator iter3=(*iter)->fillAsMuchAsPossibleWith(pol2,iter2,pol1Zip.end(),directionInPol2);
+      // and now add a full connected piece from pol1Zip:
+      if(iter3!=pol1Zip.end())
         {
           (*iter)->pushBack(*iter3);
           SoftDelete(*iter3);
-          pol2Zip.erase(iter3);
+          pol1Zip.erase(iter3);
         }
     }
 }
 
 /*!
- * 'this' is expected to be set of edges (not closed) of pol2 split.
+ * 'this' is expected to be set of edges (not closed) of pol1 split.
  */
-bool QuadraticPolygon::amIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted,const QuadraticPolygon& pol2NotSplitted, bool& direction)
+bool QuadraticPolygon::haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1NotSplitted, const QuadraticPolygon& pol2Splitted,
+                                                   bool& direction, bool& needCleaning) const
 {
-  IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(&pol1Splitted));
+  needCleaning = false;
+  IteratorOnComposedEdge it2(const_cast<QuadraticPolygon *>(&pol2Splitted));
   bool found=false;
   Node *n=getEndNode();
-  ElementaryEdge *cur=it.current();
-  for(it.first();!it.finished() && !found;)
+  ElementaryEdge *cur=it2.current();
+  // Find edge in pol2 whose start node is the end node of the current piece in pol1Zip (*this)
+  for(it2.first();!it2.finished() && !found;)
     {
-      cur=it.current();
+      cur=it2.current();
       found=(cur->getStartNode()==n);
       if(!found)
-        it.next();
+        it2.next();
     }
   if(!found)
-    throw Exception("Internal error : polygons uncompatible each others. Should never happend");
-  //Ok we found correspondance between this and pol1. Searching for right direction to close polygon.
+    throw Exception("Internal error: polygons incompatible with each others. Should never happen!");
+  //Ok we found correspondence between this and pol2. Searching for right direction to close polygon.
   ElementaryEdge *e=_sub_edges.back();
   if(e->getLoc()==FULL_ON_1)
     {
       if(e->getPtr()==cur->getPtr())
         {
-          direction=false;
-          it.previousLoop();
-          cur=it.current();
+          // if we have the same edge, several possibilities:
+          //    - either this means that pol1 and pol2 have opposite orientation (since we matched end node with start node before)
+          //    - or (more tricky, see testIntersect2DMeshes11()) that pol1 and pol2 have same orientation but 'this' turns in such
+          //    a way that it attaches to pol2 on an edge in opposite orientation.
+          // To sort this out, inspect localisation of next edge in pol2 wrt pol1NotSplitted.
+          it2.nextLoop();
+          cur=it2.current();
           Node *repr=cur->getPtr()->buildRepresentantOfMySelf();
-          bool ret=pol2NotSplitted.isInOrOut(repr);
+          bool ret=pol1NotSplitted.isInOrOut(repr);
           repr->decrRef();
+          direction = ret;
+          needCleaning = ret; // if true we are in tricky case 2 above, we know that we will produce two consecutive overlapping edges in result
           return ret;
         }
-      else
+      else  // here we don't need to go prev or next:
         {
-          direction=true;
           Node *repr=cur->getPtr()->buildRepresentantOfMySelf();
-          bool ret=pol2NotSplitted.isInOrOut(repr);
+          bool ret=pol1NotSplitted.isInOrOut(repr);
           repr->decrRef();
+          direction = ret;
           return ret;
         }
     }
   else
-    direction=cur->locateFullyMySelfAbsolute(pol2NotSplitted)==FULL_IN_1;
+    direction=cur->locateFullyMySelfAbsolute(pol1NotSplitted)==FULL_IN_1;
   return true;
 }
 
 /*!
- * This method fills as much as possible 'this' (part of pol2 split) with edges of 'pol1Splitted'.
+ * This method fills as much as possible \a this (a sub-part of pol1 split) with edges of \a pol2Splitted.
  */
-std::list<QuadraticPolygon *>::iterator QuadraticPolygon::fillAsMuchAsPossibleWith(const QuadraticPolygon& pol1Splitted,
+std::list<QuadraticPolygon *>::iterator QuadraticPolygon::fillAsMuchAsPossibleWith(const QuadraticPolygon& pol2Splitted,
                                                                                    std::list<QuadraticPolygon *>::iterator iStart,
                                                                                    std::list<QuadraticPolygon *>::iterator iEnd,
                                                                                    bool direction)
 {
-  IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(&pol1Splitted));
+  IteratorOnComposedEdge it1(const_cast<QuadraticPolygon *>(&pol2Splitted));
   bool found=false;
   Node *n=getEndNode();
-  ElementaryEdge *cur;
-  for(it.first();!it.finished() && !found;)
+  ElementaryEdge *cur1;
+  for(it1.first();!it1.finished() && !found;)
     {
-      cur=it.current();
-      found=(cur->getStartNode()==n);
+      cur1=it1.current();
+      found=(cur1->getStartNode()==n);
       if(!found)
-        it.next();
+        it1.next();
     }
   if(!direction)
-    it.previousLoop();
+    it1.previousLoop();
   Node *nodeToTest;
-  int szMax(pol1Splitted.size()+1),ii(0);// here a protection against agressive users of IntersectMeshes of invalid input meshes
+  int szMax(pol2Splitted.size()+1),ii(0);   // protection against aggressive users of IntersectMeshes using invalid input meshes ...
   std::list<QuadraticPolygon *>::iterator ret;
   do
-    {
-      cur=it.current();
-      ElementaryEdge *tmp=cur->clone();
+    { // Stack (consecutive) edges of pol1 into the result (no need to care about ordering, edges from pol1 are already consecutive)
+      cur1=it1.current();
+      ElementaryEdge *tmp=cur1->clone();
       if(!direction)
         tmp->reverse();
       pushBack(tmp);
       nodeToTest=tmp->getEndNode();
-      direction?it.nextLoop():it.previousLoop();
+      direction?it1.nextLoop():it1.previousLoop();
       ret=CheckInList(nodeToTest,iStart,iEnd);
       if(completed())
         return iEnd;
       ii++;
     }
   while(ret==iEnd && ii<szMax);
-  if(ii==szMax)// here a protection against agressive users of IntersectMeshes of invalid input meshes
+  if(ii==szMax)// here a protection against aggressive users of IntersectMeshes of invalid input meshes
     throw INTERP_KERNEL::Exception("QuadraticPolygon::fillAsMuchAsPossibleWith : Something is invalid with input polygons !");
   return ret;
 }
@@ -1151,106 +1234,157 @@ std::list<QuadraticPolygon *>::iterator QuadraticPolygon::CheckInList(Node *n, s
   return iEnd;
 }
 
-void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std::set<Edge *>& notUsedInPol1, const std::set<Edge *>& edgesInPol2OnBoundary, const std::map<INTERP_KERNEL::Node *,int>& mapp, int offset, int idThis,
-                                       std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI, std::vector<int>& nb1, std::vector<int>& nb2)
+/*!
+* Compute the remaining parts of the intersection of mesh1 by mesh2.
+* The general algorithm is to :
+*   - either return full cells from pol1 that were simply not touched by mesh2
+*   - or to:
+*      - concatenate pieces from pol1 into consecutive pieces (a bit like zipConsecutiveSegments())
+*      - complete those pieces by edges found in edgesInPol2OnBoundary, which are edges from pol2 located on the boundary of the previously built
+*      intersecting cells
+*/
+void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std::set<Edge *>& notUsedInPol1, const std::set<Edge *>& edgesInPol2OnBoundary,
+                                       const std::map<INTERP_KERNEL::Node *,mcIdType>& mapp, mcIdType offset, mcIdType idThis,
+                                       std::vector<double>& addCoordsQuadratic, std::vector<mcIdType>& conn,
+                                       std::vector<mcIdType>& connI, std::vector<mcIdType>& nb1, std::vector<mcIdType>& nb2)
 {
+  // Initialise locations on pol1. Remember that edges found in 'notUsedInPol1' are also part of the edges forming pol1.
   pol1.initLocations();
-  for(std::set<Edge *>::const_iterator it9=notUsedInPol1.begin();it9!=notUsedInPol1.end();it9++)
-    { (*it9)->initLocs(); (*it9)->declareOn(); }
-  for(std::set<Edge *>::const_iterator itA=edgesInPol2OnBoundary.begin();itA!=edgesInPol2OnBoundary.end();itA++)
-    { (*itA)->initLocs(); (*itA)->declareIn(); }
+  for(std::set<Edge *>::const_iterator it1=notUsedInPol1.begin();it1!=notUsedInPol1.end();it1++)
+    { (*it1)->initLocs(); (*it1)->declareOn(); }
+  for(std::set<Edge *>::const_iterator it2=edgesInPol2OnBoundary.begin();it2!=edgesInPol2OnBoundary.end();it2++)
+    { (*it2)->initLocs(); (*it2)->declareIn(); }
   ////
   std::set<Edge *> notUsedInPol1L(notUsedInPol1);
-  IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(&pol1));
+  IteratorOnComposedEdge itPol1(const_cast<QuadraticPolygon *>(&pol1));
   int sz=pol1.size();
   std::list<QuadraticPolygon *> pol1Zip;
+  // If none of the edges of pol1 was consumed by the rebuilding process, we can directly take pol1 as it is to form a cell:
   if(pol1.size()==(int)notUsedInPol1.size() && edgesInPol2OnBoundary.empty())
     {
       pol1.appendCrudeData(mapp,0.,0.,1.,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1);
       return ;
     }
+  // Zip consecutive edges found in notUsedInPol1L and which are not overlapping with boundary edge from edgesInPol2OnBoundary:
   while(!notUsedInPol1L.empty())
     {
-      for(int i=0;i<sz && (it.current()->getStartNode()->getLoc()!=IN_1 || it.current()->getLoc()!=FULL_ON_1);i++)
-        it.nextLoop();
-      if(it.current()->getStartNode()->getLoc()!=IN_1 || it.current()->getLoc()!=FULL_ON_1)
+      // If all nodes are IN or edge is ON (i.e. as initialised at the begining of the method) then error
+      for(int i=0;i<sz && (itPol1.current()->getStartNode()->getLoc()!=IN_1 || itPol1.current()->getLoc()!=FULL_ON_1);i++)
+        itPol1.nextLoop();
+      if(itPol1.current()->getStartNode()->getLoc()!=IN_1 || itPol1.current()->getLoc()!=FULL_ON_1)
         throw INTERP_KERNEL::Exception("Presence of a target polygon fully included in source polygon ! The partition of this leads to a non simply connex cell (with hole) ! Impossible ! Such resulting cell cannot be stored in MED cell format !");
       QuadraticPolygon *tmp1=new QuadraticPolygon;
       do
         {
-          Edge *ee=it.current()->getPtr();
+          Edge *ee=itPol1.current()->getPtr();
           if(ee->getLoc()==FULL_ON_1)
             {
               ee->incrRef(); notUsedInPol1L.erase(ee);
-              tmp1->pushBack(new ElementaryEdge(ee,it.current()->getDirection()));    
+              tmp1->pushBack(new ElementaryEdge(ee,itPol1.current()->getDirection()));
             }
-          it.nextLoop();
+          itPol1.nextLoop();
         }
-      while(it.current()->getStartNode()->getLoc()!=IN_1 && !notUsedInPol1L.empty());
+      while(itPol1.current()->getStartNode()->getLoc()!=IN_1 && !notUsedInPol1L.empty());
       pol1Zip.push_back(tmp1);
     }
+
   ////
   std::list<QuadraticPolygon *> retPolsUnderContruction;
   std::list<Edge *> edgesInPol2OnBoundaryL(edgesInPol2OnBoundary.begin(),edgesInPol2OnBoundary.end());
-  std::map<QuadraticPolygon *, std::list<QuadraticPolygon *> > pol1ZipConsumed;
+  std::map<QuadraticPolygon *, std::list<QuadraticPolygon *> > pol1ZipConsumed;  // for memory management only.
   std::size_t maxNbOfTurn=edgesInPol2OnBoundaryL.size(),nbOfTurn=0,iiMNT=0;
   for(std::list<QuadraticPolygon *>::const_iterator itMNT=pol1Zip.begin();itMNT!=pol1Zip.end();itMNT++,iiMNT++)
     nbOfTurn+=(*itMNT)->size();
   maxNbOfTurn=maxNbOfTurn*nbOfTurn; maxNbOfTurn*=maxNbOfTurn;
+  // [ABN] at least 3 turns for very small cases (typically one (quad) edge against one (quad or lin) edge forming a new cell)!
+  maxNbOfTurn = maxNbOfTurn<3 ? 3 : maxNbOfTurn;
   nbOfTurn=0;
-  while(nbOfTurn<maxNbOfTurn && ((!pol1Zip.empty() || !edgesInPol2OnBoundaryL.empty())))
+  while(nbOfTurn<maxNbOfTurn)  // the 'normal' way out of this loop is the break towards the end when pol1Zip is empty.
     {
-      for(std::list<QuadraticPolygon *>::iterator it1=retPolsUnderContruction.begin();it1!=retPolsUnderContruction.end();)
+      // retPolsUnderConstruction initially empty -> see if(!pol1Zip.empty()) below ...
+      for(std::list<QuadraticPolygon *>::iterator itConstr=retPolsUnderContruction.begin();itConstr!=retPolsUnderContruction.end();)
         {
-          if((*it1)->getStartNode()==(*it1)->getEndNode())
-            {
-              it1++;
-              continue;
-            }
-          Node *curN=(*it1)->getEndNode();
-          bool smthHappened=false;
+          Node *startN = (*itConstr)->getStartNode();
+          Node *curN = (*itConstr)->getEndNode();
+          if(startN == curN)  // reconstruction of a cell is finished
+            {  itConstr++; continue;  }
+
+          bool smthHappened=false, doneEarly=false;
+          // Complete a partially reconstructed polygon with boundary edges of pol2 by matching nodes:
           for(std::list<Edge *>::iterator it2=edgesInPol2OnBoundaryL.begin();it2!=edgesInPol2OnBoundaryL.end();)
             {
-              if(curN==(*it2)->getStartNode())
-                { (*it2)->incrRef(); (*it1)->pushBack(new ElementaryEdge(*it2,true)); curN=(*it2)->getEndNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); }
-              else if(curN==(*it2)->getEndNode())
-                { (*it2)->incrRef(); (*it1)->pushBack(new ElementaryEdge(*it2,false)); curN=(*it2)->getStartNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); }
+              if(curN==(*it2)->getEndNode())  // only end node should be considered if orientation is correct for input meshes
+                                              // in the funny case of cells exactly included (see test case testIntersect2DMeshesTmp13) this is mandatory to take edges from pol2 in the right order.
+                {
+                  (*it2)->incrRef();
+                  (*itConstr)->pushBack(new ElementaryEdge(*it2,false));
+                  curN=(*it2)->getStartNode();
+                  smthHappened=true;
+                  it2=edgesInPol2OnBoundaryL.erase(it2);
+                }
               else
                 it2++;
+              if (curN == startN) // we might end here
+                { doneEarly = true;  break;  }
             }
-          if(smthHappened)
+
+          // It might be the case that the lookup on start nodes made above failed because pol2 is wrongly oriented.
+          // Be somewhat flexible and keep on supporting this case here (useful for voronisation notably):
+          if(!smthHappened)
             {
-              for(std::list<QuadraticPolygon *>::iterator it3=pol1Zip.begin();it3!=pol1Zip.end();)
+              for(std::list<Edge *>::iterator it2=edgesInPol2OnBoundaryL.begin();it2!=edgesInPol2OnBoundaryL.end();)
                 {
-                  if(curN==(*it3)->getStartNode())
+                  if(curN==(*it2)->getStartNode())
                     {
-                      for(std::list<ElementaryEdge *>::const_iterator it4=(*it3)->_sub_edges.begin();it4!=(*it3)->_sub_edges.end();it4++)
-                        { (*it4)->getPtr()->incrRef(); bool dir=(*it4)->getDirection(); (*it1)->pushBack(new ElementaryEdge((*it4)->getPtr(),dir)); }
+                      (*it2)->incrRef();
+                      (*itConstr)->pushBack(new ElementaryEdge(*it2,true));
+                      curN=(*it2)->getEndNode();
                       smthHappened=true;
-                      pol1ZipConsumed[*it1].push_back(*it3);
-                      curN=(*it3)->getEndNode();
-                      it3=pol1Zip.erase(it3);
+                      it2=edgesInPol2OnBoundaryL.erase(it2);
                     }
                   else
-                    it3++;
+                    it2++;
+                  if (curN == startN) // we might end here
+                      { doneEarly = true;  break;  }
                 }
             }
-          if(!smthHappened)
+
+          if (doneEarly)
+            {  itConstr++; continue;  }
+
+          if(smthHappened) // Now continue the construction by finding the next bit in pol1Zip. Not too sure what are the cases where the boolean if False ...
             {
-              for(std::list<ElementaryEdge *>::const_iterator it5=(*it1)->_sub_edges.begin();it5!=(*it1)->_sub_edges.end();it5++)
+              for(std::list<QuadraticPolygon *>::iterator itZip=pol1Zip.begin();itZip!=pol1Zip.end();)
+                {
+                  if(curN==(*itZip)->getStartNode()) // we found a matching piece to append in pol1Zip. Append all of it to the current polygon being reconstr
+                    {
+                      for(std::list<ElementaryEdge *>::const_iterator it4=(*itZip)->_sub_edges.begin();it4!=(*itZip)->_sub_edges.end();it4++)
+                        { (*it4)->getPtr()->incrRef(); bool dir=(*it4)->getDirection(); (*itConstr)->pushBack(new ElementaryEdge((*it4)->getPtr(),dir)); }
+                      pol1ZipConsumed[*itConstr].push_back(*itZip);
+                      curN=(*itZip)->getEndNode();
+                      itZip=pol1Zip.erase(itZip);  // one zipped piece has been consumed
+                      break;                       // we can stop here, pieces in pol1Zip are not connected, by definition.
+                    }
+                  else
+                    itZip++;
+                }
+            }
+          else // Nothing happened.
+            {
+              for(std::list<ElementaryEdge *>::const_iterator it5=(*itConstr)->_sub_edges.begin();it5!=(*itConstr)->_sub_edges.end();it5++)
                 {
                   Edge *ee=(*it5)->getPtr();
                   if(edgesInPol2OnBoundary.find(ee)!=edgesInPol2OnBoundary.end())
                     edgesInPol2OnBoundaryL.push_back(ee);
                 }
-              for(std::list<QuadraticPolygon *>::iterator it6=pol1ZipConsumed[*it1].begin();it6!=pol1ZipConsumed[*it1].end();it6++)
+              for(std::list<QuadraticPolygon *>::iterator it6=pol1ZipConsumed[*itConstr].begin();it6!=pol1ZipConsumed[*itConstr].end();it6++)
                 pol1Zip.push_front(*it6);
-              pol1ZipConsumed.erase(*it1);
-              delete *it1;
-              it1=retPolsUnderContruction.erase(it1);
+              pol1ZipConsumed.erase(*itConstr);
+              delete *itConstr;
+              itConstr=retPolsUnderContruction.erase(itConstr);
             }
         }
-      if(!pol1Zip.empty())
+      if(!pol1Zip.empty())  // the filling process of retPolsUnderConstruction starts here
         {
           QuadraticPolygon *tmp=new QuadraticPolygon;
           QuadraticPolygon *first=*(pol1Zip.begin());
@@ -1260,6 +1394,8 @@ void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std::
           retPolsUnderContruction.push_back(tmp);
           pol1Zip.erase(pol1Zip.begin());
         }
+      else
+        break;
       nbOfTurn++;
     }
   if(nbOfTurn==maxNbOfTurn)
@@ -1268,14 +1404,21 @@ void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std::
       oss << " Number of turns is = " << nbOfTurn << " !";
       throw INTERP_KERNEL::Exception(oss.str().c_str());
     }
-  for(std::list<QuadraticPolygon *>::iterator it1=retPolsUnderContruction.begin();it1!=retPolsUnderContruction.end();it1++)
+  // Convert to integer connectivity:
+  for(std::list<QuadraticPolygon *>::iterator itConstr=retPolsUnderContruction.begin();itConstr!=retPolsUnderContruction.end();itConstr++)
     {
-      if((*it1)->getStartNode()==(*it1)->getEndNode())
+      if((*itConstr)->getStartNode()==(*itConstr)->getEndNode())  // take only fully closed reconstructed polygon
         {
-          (*it1)->appendCrudeData(mapp,0.,0.,1.,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1);
-          for(std::list<QuadraticPolygon *>::iterator it6=pol1ZipConsumed[*it1].begin();it6!=pol1ZipConsumed[*it1].end();it6++)
+          (*itConstr)->cleanDegeneratedConsecutiveEdges();
+          (*itConstr)->appendCrudeData(mapp,0.,0.,1.,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1);
+          for(std::list<QuadraticPolygon *>::iterator it6=pol1ZipConsumed[*itConstr].begin();it6!=pol1ZipConsumed[*itConstr].end();it6++)
             delete *it6;
-          delete *it1;
+          delete *itConstr;
+        }
+      else
+        {
+          std::ostringstream oss; oss << "Internal error during reconstruction of residual of cell! Non fully closed polygon built!";
+          throw INTERP_KERNEL::Exception(oss.str().c_str());
         }
     }
 }