src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
  22 #include "InterpKernelGeo2DElementaryEdge.hxx"
  23 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
  24 #include "InterpKernelGeo2DAbstractEdge.hxx"
  25 #include "InterpKernelGeo2DEdgeLin.hxx"
  26 #include "InterpKernelGeo2DBounds.hxx"
  27 #include "InterpKernelGeo2DEdge.txx"
  28
  29 #include "NormalizedUnstructuredMesh.hxx"
  30
  31 #include <fstream>
  32 #include <sstream>
  33 #include <iomanip>
  34 #include <cstring>
  35 #include <limits>
  36
  37 using namespace INTERP_KERNEL;
  38
  39 namespace INTERP_KERNEL
  40 {
  41   const unsigned MAX_SIZE_OF_LINE_XFIG_FILE=1024;
  42 }
  43
  44 QuadraticPolygon::QuadraticPolygon(const char *file)
  45 {
  46   char currentLine[MAX_SIZE_OF_LINE_XFIG_FILE];
  47   std::ifstream stream(file);
  48   stream.exceptions(std::ios_base::eofbit);
  49   try
  50     {
  51       do
  52         stream.getline(currentLine,MAX_SIZE_OF_LINE_XFIG_FILE);
  53       while(strcmp(currentLine,"1200 2")!=0);
  54       do
  55         {
  56           Edge *newEdge=Edge::BuildFromXfigLine(stream);
  57           if(!empty())
  58             newEdge->changeStartNodeWith(back()->getEndNode());
  59           pushBack(newEdge);
  60         }
  61       while(1);
  62     }
  63   catch(std::ifstream::failure&)
  64     {
  65     }
  66   front()->changeStartNodeWith(back()->getEndNode());
  67 }
  68
  69 QuadraticPolygon::~QuadraticPolygon()
  70 {
  71 }
  72
  73 QuadraticPolygon *QuadraticPolygon::BuildLinearPolygon(std::vector<Node *>& nodes, bool isClosed)
  74 {
  75   QuadraticPolygon *ret=new QuadraticPolygon;
  76   std::size_t size=nodes.size();
  77   for(std::size_t i=0;i< (size - (isClosed?0:1));i++)
  78     {
  79       ret->pushBack(new EdgeLin(nodes[i],nodes[(i+1)%size]));
  80       nodes[i]->decrRef();
  81     }
  82   if(!isClosed)
  83     nodes[size-1]->decrRef();
  84   return ret;
  85 }
  86
  87 QuadraticPolygon *QuadraticPolygon::BuildArcCirclePolygon(std::vector<Node *>& nodes, bool isClosed)
  88 {
  89   QuadraticPolygon *ret=new QuadraticPolygon;
  90   std::size_t size=nodes.size();
  91   std::size_t quad_offset = isClosed ? (size/2) : (size/2+1);
  92   for(std::size_t i = 0; i < size/2; i++)
  93     {
  94       EdgeLin *e1,*e2;
  95       e1=new EdgeLin(nodes[i],nodes[i+quad_offset]);
  96       e2=new EdgeLin(nodes[i+quad_offset],nodes[(i+1)%quad_offset]);
  97       SegSegIntersector inters(*e1,*e2);
  98       bool colinearity=inters.areColinears();
  99       delete e1; delete e2;
 100       if(colinearity)
 101         ret->pushBack(new EdgeLin(nodes[i],nodes[(i+1)%quad_offset]));
 102       else
 103         ret->pushBack(new EdgeArcCircle(nodes[i],nodes[i+quad_offset],nodes[(i+1)%quad_offset]));
 104       nodes[i]->decrRef(); nodes[i+quad_offset]->decrRef();
 105     }
 106   if(!isClosed)
 107     nodes[quad_offset-1]->decrRef();
 108   return ret;
 109 }
 110
 111 void QuadraticPolygon::BuildDbgFile(const std::vector<Node *>& nodes, const char *fileName)
 112 {
 113   std::ofstream file(fileName);
 114   file << std::setprecision(16);
 115   file << "  double coords[]=" << std::endl << "    { ";
 116   for(std::vector<Node *>::const_iterator iter=nodes.begin();iter!=nodes.end();iter++)
 117     {
 118       if(iter!=nodes.begin())
 119         file << "," << std::endl << "      ";
 120       file << (*(*iter))[0] << ", " << (*(*iter))[1];
 121     }
 122   file << "};" << std::endl;
 123 }
 124
 125 void QuadraticPolygon::closeMe() const
 126 {
 127   if(!front()->changeStartNodeWith(back()->getEndNode()))
 128     throw(Exception("big error: not closed polygon..."));
 129 }
 130
 131 void QuadraticPolygon::circularPermute()
 132 {
 133   if(_sub_edges.size()>1)
 134     {
 135       ElementaryEdge *first=_sub_edges.front();
 136       _sub_edges.pop_front();
 137       _sub_edges.push_back(first);
 138     }
 139 }
 140
 141 bool QuadraticPolygon::isButterflyAbs()
 142 {
 143   INTERP_KERNEL::Bounds b;
 144   double xBary,yBary;
 145   b.prepareForAggregation();
 146   fillBounds(b);
 147   double dimChar=b.getCaracteristicDim();
 148   b.getBarycenter(xBary,yBary);
 149   applyGlobalSimilarity(xBary,yBary,dimChar);
 150   //
 151   return isButterfly();
 152 }
 153
 154 bool QuadraticPolygon::isButterfly() const
 155 {
 156   for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++)
 157     {
 158       Edge *e1=(*it)->getPtr();
 159       std::list<ElementaryEdge *>::const_iterator it2=it;
 160       it2++;
 161       for(;it2!=_sub_edges.end();it2++)
 162         {
 163           MergePoints commonNode;
 164           ComposedEdge *outVal1=new ComposedEdge;
 165           ComposedEdge *outVal2=new ComposedEdge;
 166           Edge *e2=(*it2)->getPtr();
 167           if(e1->intersectWith(e2,commonNode,*outVal1,*outVal2))
 168             {
 169               Delete(outVal1);
 170               Delete(outVal2);
 171               return true;
 172             }
 173           Delete(outVal1);
 174           Delete(outVal2);
 175         }
 176     }
 177   return false;
 178 }
 179
 180 void QuadraticPolygon::dumpInXfigFileWithOther(const ComposedEdge& other, const char *fileName) const
 181 {
 182   std::ofstream file(fileName);
 183   const int resolution=1200;
 184   Bounds box;
 185   box.prepareForAggregation();
 186   fillBounds(box);
 187   other.fillBounds(box);
 188   dumpInXfigFile(file,resolution,box);
 189   other.ComposedEdge::dumpInXfigFile(file,resolution,box);
 190 }
 191
 192 void QuadraticPolygon::dumpInXfigFile(const char *fileName) const
 193 {
 194   std::ofstream file(fileName);
 195   const int resolution=1200;
 196   Bounds box;
 197   box.prepareForAggregation();
 198   fillBounds(box);
 199   dumpInXfigFile(file,resolution,box);
 200 }
 201
 202 void QuadraticPolygon::dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const
 203 {
 204   stream << "#FIG 3.2  Produced by xfig version 3.2.5-alpha5" << std::endl;
 205   stream << "Landscape" << std::endl;
 206   stream << "Center" << std::endl;
 207   stream << "Metric" << std::endl;
 208   stream << "Letter" << std::endl;
 209   stream << "100.00" << std::endl;
 210   stream << "Single" << std::endl;
 211   stream << "-2" << std::endl;
 212   stream << resolution << " 2" << std::endl;
 213   ComposedEdge::dumpInXfigFile(stream,resolution,box);
 214 }
 215
 216 /*!
 217  * Warning contrary to intersectWith method this method is \b NOT const. 'this' and 'other' are modified after call of this method.
 218  */
 219 double QuadraticPolygon::intersectWithAbs(QuadraticPolygon& other)
 220 {
 221   double ret=0.,xBaryBB,yBaryBB;
 222   double fact=normalize(&other,xBaryBB,yBaryBB);
 223   std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
 224   for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
 225     {
 226       ret+=fabs((*iter)->getArea());
 227       delete *iter;
 228     }
 229   return ret*fact*fact;
 230 }
 231
 232 /*!
 233  * This method splits 'this' with 'other' into smaller pieces localizable. 'mapThis' is a map that gives the correspondance
 234  * between nodes contained in 'this' and node ids in a global mesh.
 235  * In the same way, 'mapOther' gives the correspondance between nodes contained in 'other' and node ids in a
 236  * global mesh from wich 'other' is extracted.
 237  * This method has 1 out paramater : 'edgesThis', After the call of this method, it contains the nodal connectivity (including type)
 238  * of 'this' into globlal "this mesh".
 239  * This method has 2 in/out parameters : 'subDivOther' and 'addCoo'.'otherEdgeIds' is useful to put values in
 240  * 'edgesThis', 'subDivOther' and 'addCoo'.
 241  * Size of 'otherEdgeIds' has to be equal to number of ElementaryEdges in 'other'. No check of that will be done.
 242  * The term 'abs' in the name recalls that we normalize the mesh (spatially) so that node coordinates fit into [0;1].
 243  * @param offset1 is the number of nodes contained in global mesh from which 'this' is extracted.
 244  * @param offset2 is the sum of nodes contained in global mesh from which 'this' is extracted and 'other' is extracted.
 245  * @param edgesInOtherColinearWithThis will be appended at the end of the vector with colinear edge ids of other (if any)
 246  * @param otherEdgeIds is a vector with the same size than other before calling this method. It gives in the same order
 247  * the cell id in global other mesh.
 248  */
 249 void QuadraticPolygon::splitAbs(QuadraticPolygon& other,
 250         const std::map<INTERP_KERNEL::Node *,int>& mapThis, const std::map<INTERP_KERNEL::Node *,int>& mapOther,
 251         int offset1, int offset2 ,
 252         const std::vector<int>& otherEdgeIds,
 253         std::vector<int>& edgesThis, int cellIdThis,
 254         std::vector< std::vector<int> >& edgesInOtherColinearWithThis, std::vector< std::vector<int> >& subDivOther,
 255         std::vector<double>& addCoo)
 256 {
 257   double xBaryBB, yBaryBB;
 258   double fact=normalizeExt(&other, xBaryBB, yBaryBB);
 259   //
 260   IteratorOnComposedEdge it1(this),it3(&other);
 261   MergePoints merge;
 262   ComposedEdge *c1=new ComposedEdge;
 263   ComposedEdge *c2=new ComposedEdge;
 264   int i=0;
 265   std::map<INTERP_KERNEL::Node *,int> mapAddCoo;
 266   for(it3.first();!it3.finished();it3.next(),i++)//iteration over 'other' _sub_edges
 267     {
 268       QuadraticPolygon otherTmp;
 269       ElementaryEdge* curE3=it3.current();
 270       otherTmp.pushBack(new ElementaryEdge(curE3->getPtr(),curE3->getDirection())); curE3->getPtr()->incrRef();
 271       IteratorOnComposedEdge it2(&otherTmp);
 272       for(it2.first();!it2.finished();it2.next())//iteration on subedges of 'other->_sub_edge'
 273         {
 274           ElementaryEdge* curE2=it2.current();
 275           if(!curE2->isThereStartPoint())
 276             it1.first();
 277           else
 278             it1=curE2->getIterator();
 279           for(;!it1.finished();)//iteration over 'this' _sub_edges
 280             {
 281               ElementaryEdge* curE1=it1.current();
 282               merge.clear();
 283               if(curE1->getPtr()->intersectWith(curE2->getPtr(),merge,*c1,*c2))
 284                 {
 285                   if(!curE1->getDirection()) c1->reverse();
 286                   if(!curE2->getDirection()) c2->reverse();
 287                   UpdateNeighbours(merge,it1,it2,c1,c2);
 288                   //Substitution of simple edge by sub-edges.
 289                   delete curE1; // <-- destroying simple edge coming from pol1
 290                   delete curE2; // <-- destroying simple edge coming from pol2
 291                   it1.insertElemEdges(c1,true);// <-- 2nd param is true to go next.
 292                   it2.insertElemEdges(c2,false);// <-- 2nd param is false to avoid to go next.
 293                   curE2=it2.current();
 294                   //
 295                   it1.assignMySelfToAllElems(c2);//To avoid that others
 296                   SoftDelete(c1);
 297                   SoftDelete(c2);
 298                   c1=new ComposedEdge;
 299                   c2=new ComposedEdge;
 300                 }
 301               else
 302                 {
 303                   UpdateNeighbours(merge,it1,it2,curE1,curE2);
 304                   it1.next();
 305                 }
 306             }
 307         }
 308       if(otherTmp.presenceOfOn())
 309         edgesInOtherColinearWithThis[otherEdgeIds[i]].push_back(cellIdThis);
 310       if(otherTmp._sub_edges.size()>1)
 311         {
 312           for(std::list<ElementaryEdge *>::const_iterator it=otherTmp._sub_edges.begin();it!=otherTmp._sub_edges.end();it++)
 313             (*it)->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,/**/fact,xBaryBB,yBaryBB,/**/subDivOther[otherEdgeIds[i]],addCoo,mapAddCoo);
 314         }
 315     }
 316   Delete(c1);
 317   Delete(c2);
 318   //
 319   for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++)
 320     (*it)->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,/**/fact,xBaryBB,yBaryBB,/**/edgesThis,addCoo,mapAddCoo);
 321   //
 322 }
 323
 324 /*!
 325  * This method builds 'this' from its descending conn stored in crude mode (MEDCoupling).
 326  * Descending conn is in FORTRAN relative mode in order to give the
 327  * orientation of edge (see buildDescendingConnectivity2() method).
 328  * See appendEdgeFromCrudeDataArray() for params description.
 329  */
 330 void QuadraticPolygon::buildFromCrudeDataArray(const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad, const int *nodalBg, const double *coords,
 331                                                const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges)
 332 {
 333   std::size_t nbOfSeg=std::distance(descBg,descEnd);
 334   for(std::size_t i=0;i<nbOfSeg;i++)
 335     {
 336       appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
 337     }
 338 }
 339
 340 void QuadraticPolygon::appendEdgeFromCrudeDataArray(std::size_t edgePos, const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad,
 341                             const int *nodalBg, const double *coords,
 342                             const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges)
 343 {
 344   if(!isQuad)
 345     {
 346       bool direct=descBg[edgePos]>0;
 347       int edgeId=abs(descBg[edgePos])-1; // back to C indexing mode
 348       const std::vector<int>& subEdge=intersectEdges[edgeId];
 349       std::size_t nbOfSubEdges=subEdge.size()/2;
 350       for(std::size_t j=0;j<nbOfSubEdges;j++)
 351         appendSubEdgeFromCrudeDataArray(0,j,direct,edgeId,subEdge,mapp);
 352     }
 353   else
 354     {
 355       std::size_t nbOfSeg=std::distance(descBg,descEnd);
 356       const double *st=coords+2*(nodalBg[edgePos]);
 357       INTERP_KERNEL::Node *st0=new INTERP_KERNEL::Node(st[0],st[1]);
 358       const double *endd=coords+2*(nodalBg[(edgePos+1)%nbOfSeg]);
 359       INTERP_KERNEL::Node *endd0=new INTERP_KERNEL::Node(endd[0],endd[1]);
 360       const double *middle=coords+2*(nodalBg[edgePos+nbOfSeg]);
 361       INTERP_KERNEL::Node *middle0=new INTERP_KERNEL::Node(middle[0],middle[1]);
 362       EdgeLin *e1,*e2;
 363       e1=new EdgeLin(st0,middle0);
 364       e2=new EdgeLin(middle0,endd0);
 365       SegSegIntersector inters(*e1,*e2);
 366       bool colinearity=inters.areColinears();
 367       delete e1; delete e2;
 368       //
 369       bool direct=descBg[edgePos]>0;
 370       int edgeId=abs(descBg[edgePos])-1;
 371       const std::vector<int>& subEdge=intersectEdges[edgeId];
 372       std::size_t nbOfSubEdges=subEdge.size()/2;
 373       if(colinearity)
 374         {
 375           for(std::size_t j=0;j<nbOfSubEdges;j++)
 376             appendSubEdgeFromCrudeDataArray(0,j,direct,edgeId,subEdge,mapp);
 377         }
 378       else
 379         {
 380           Edge *e=new EdgeArcCircle(st0,middle0,endd0,true);
 381           for(std::size_t j=0;j<nbOfSubEdges;j++)
 382             appendSubEdgeFromCrudeDataArray(e,j,direct,edgeId,subEdge,mapp);
 383           e->decrRef();
 384         }
 385       st0->decrRef(); endd0->decrRef(); middle0->decrRef();
 386     }
 387 }
 388
 389 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)
 390 {
 391   std::size_t nbOfSubEdges=subEdge.size()/2;
 392   if(!baseEdge)
 393     {//it is not a quadratic subedge
 394       Node *start=(*mapp.find(direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1])).second;
 395       Node *end=(*mapp.find(direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2])).second;
 396       ElementaryEdge *e=ElementaryEdge::BuildEdgeFromCrudeDataArray(true,start,end);
 397       pushBack(e);
 398     }
 399   else
 400     {//it is a quadratic subedge
 401       Node *start=(*mapp.find(direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1])).second;
 402       Node *end=(*mapp.find(direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2])).second;
 403       Edge *ee=baseEdge->buildEdgeLyingOnMe(start,end);
 404       ElementaryEdge *eee=new ElementaryEdge(ee,true);
 405       pushBack(eee);
 406     }
 407 }
 408
 409 /*!
 410  * 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
 411  * orientation of edge.
 412  */
 413 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,
 414                                                 const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, const std::vector<std::vector<int> >& intersectEdges1,
 415                                                 const std::vector< std::vector<int> >& colinear1,
 416                                                 std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> >& alreadyExistingIn2)
 417 {
 418   std::size_t nbOfSeg=std::distance(descBg,descEnd);
 419   for(std::size_t i=0;i<nbOfSeg;i++)//loop over all edges of pol2
 420     {
 421       bool direct=descBg[i]>0;
 422       int edgeId=abs(descBg[i])-1;//current edge id of pol2
 423       std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> >::const_iterator it1=alreadyExistingIn2.find(descBg[i]),it2=alreadyExistingIn2.find(-descBg[i]);
 424       if(it1!=alreadyExistingIn2.end() || it2!=alreadyExistingIn2.end())
 425         {
 426           bool sameDir=(it1!=alreadyExistingIn2.end());
 427           const std::vector<INTERP_KERNEL::ElementaryEdge *>& edgesAlreadyBuilt=sameDir?(*it1).second:(*it2).second;
 428           if(sameDir)
 429             {
 430               for(std::vector<INTERP_KERNEL::ElementaryEdge *>::const_iterator it3=edgesAlreadyBuilt.begin();it3!=edgesAlreadyBuilt.end();it3++)
 431                 {
 432                   Edge *ee=(*it3)->getPtr(); ee->incrRef();
 433                   pushBack(new ElementaryEdge(ee,(*it3)->getDirection()));
 434                 }
 435             }
 436           else
 437             {
 438               for(std::vector<INTERP_KERNEL::ElementaryEdge *>::const_reverse_iterator it4=edgesAlreadyBuilt.rbegin();it4!=edgesAlreadyBuilt.rend();it4++)
 439                 {
 440                   Edge *ee=(*it4)->getPtr(); ee->incrRef();
 441                   pushBack(new ElementaryEdge(ee,!(*it4)->getDirection()));
 442                 }
 443             }
 444           continue;
 445         }
 446       bool directos=colinear1[edgeId].empty();
 447       std::vector<std::pair<int,std::pair<bool,int> > > idIns1;
 448       int offset1=0;
 449       if(!directos)
 450         {// if the current edge of pol2 has one or more colinear edges part into pol1
 451           const std::vector<int>& c=colinear1[edgeId];
 452           std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1);
 453           for(std::size_t j=0;j<nbOfEdgesIn1;j++)
 454             {
 455               int edgeId1=abs(descBg1[j])-1;
 456               if(std::find(c.begin(),c.end(),edgeId1)!=c.end())
 457                 {
 458                   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
 459                   //std::pair<edgeId1); direct1=descBg1[j]>0;
 460                 }
 461               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
 462             }
 463           directos=idIns1.empty();
 464         }
 465       if(directos)
 466         {//no subpart of edge 'edgeId' of pol2 is in pol1 so let's operate the same thing that QuadraticPolygon::buildFromCrudeDataArray method
 467           std::size_t oldSz=_sub_edges.size();
 468           appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
 469           std::size_t newSz=_sub_edges.size();
 470           std::size_t zeSz=newSz-oldSz;
 471           alreadyExistingIn2[descBg[i]].resize(zeSz);
 472           std::list<ElementaryEdge *>::const_reverse_iterator it5=_sub_edges.rbegin();
 473           for(std::size_t p=0;p<zeSz;p++,it5++)
 474             alreadyExistingIn2[descBg[i]][zeSz-p-1]=*it5;
 475         }
 476       else
 477         {//there is subpart of edge 'edgeId' of pol2 inside pol1
 478           const std::vector<int>& subEdge=intersectEdges[edgeId];
 479           std::size_t nbOfSubEdges=subEdge.size()/2;
 480           for(std::size_t j=0;j<nbOfSubEdges;j++)
 481             {
 482               int idBg=direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1];
 483               int idEnd=direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2];
 484               bool direction11,found=false;
 485               bool direct1;//store if needed the direction in 1
 486               int offset2;
 487               std::size_t nbOfSubEdges1;
 488               for(std::vector<std::pair<int,std::pair<bool,int> > >::const_iterator it=idIns1.begin();it!=idIns1.end() && !found;it++)
 489                 {
 490                   int idIn1=(*it).first;//store if needed the cell id in 1
 491                   direct1=(*it).second.first;
 492                   offset1=(*it).second.second;
 493                   const std::vector<int>& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1];
 494                   nbOfSubEdges1=subEdge1PossiblyAlreadyIn1.size()/2;
 495                   offset2=0;
 496                   for(std::size_t k=0;k<nbOfSubEdges1 && !found;k++)
 497                     {//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
 498                       if(subEdge1PossiblyAlreadyIn1[2*k]==idBg && subEdge1PossiblyAlreadyIn1[2*k+1]==idEnd)
 499                         { direction11=true; found=true; }
 500                       else if(subEdge1PossiblyAlreadyIn1[2*k]==idEnd && subEdge1PossiblyAlreadyIn1[2*k+1]==idBg)
 501                         { direction11=false; found=true; }
 502                       else
 503                         offset2++;
 504                     }
 505                 }
 506               if(!found)
 507                 {//the current subedge of edge 'edgeId' of pol2 is not a part of the colinear edge 'idIn1' of pol1 -> build new Edge instance
 508                   //appendEdgeFromCrudeDataArray(j,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
 509                   Node *start=(*mapp.find(idBg)).second;
 510                   Node *end=(*mapp.find(idEnd)).second;
 511                   ElementaryEdge *e=ElementaryEdge::BuildEdgeFromCrudeDataArray(true,start,end);
 512                   pushBack(e);
 513                   alreadyExistingIn2[descBg[i]].push_back(e);
 514                 }
 515               else
 516                 {//the current subedge of edge 'edgeId' of pol2 is part of the colinear edge 'idIn1' of pol1 -> reuse Edge instance of pol1
 517                   ElementaryEdge *e=pol1[offset1+(direct1?offset2:nbOfSubEdges1-offset2-1)];
 518                   Edge *ee=e->getPtr();
 519                   ee->incrRef();
 520                   ElementaryEdge *e2=new ElementaryEdge(ee,!(direct1^direction11));
 521                   pushBack(e2);
 522                   alreadyExistingIn2[descBg[i]].push_back(e2);
 523                 }
 524             }
 525         }
 526     }
 527 }
 528
 529 /*!
 530  * Method expected to be called on pol2. Every params not suffixed by numbered are supposed to refer to pol2 (this).
 531  * Method to find edges that are ON.
 532  */
 533 void QuadraticPolygon::updateLocOfEdgeFromCrudeDataArray2(const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges,
 534       const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1,
 535       const std::vector<std::vector<int> >& intersectEdges1, const std::vector< std::vector<int> >& colinear1) const
 536 {
 537   std::size_t nbOfSeg=std::distance(descBg,descEnd);
 538   for(std::size_t i=0;i<nbOfSeg;i++)//loop over all edges of pol2
 539     {
 540       bool direct=descBg[i]>0;
 541       int edgeId=abs(descBg[i])-1;//current edge id of pol2
 542       const std::vector<int>& c=colinear1[edgeId];
 543       if(c.empty())
 544         continue;
 545       const std::vector<int>& subEdge=intersectEdges[edgeId];
 546       std::size_t nbOfSubEdges=subEdge.size()/2;
 547       //
 548       std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1);
 549       int offset1=0;
 550       for(std::size_t j=0;j<nbOfEdgesIn1;j++)
 551         {
 552           int edgeId1=abs(descBg1[j])-1;
 553           if(std::find(c.begin(),c.end(),edgeId1)!=c.end())
 554             {
 555               for(std::size_t k=0;k<nbOfSubEdges;k++)
 556                 {
 557                   int idBg=direct?subEdge[2*k]:subEdge[2*nbOfSubEdges-2*k-1];
 558                   int idEnd=direct?subEdge[2*k+1]:subEdge[2*nbOfSubEdges-2*k-2];
 559                   int idIn1=edgeId1;
 560                   bool direct1=descBg1[j]>0;
 561                   const std::vector<int>& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1];
 562                   std::size_t nbOfSubEdges1=subEdge1PossiblyAlreadyIn1.size()/2;
 563                   int offset2=0;
 564                   bool found=false;
 565                   for(std::size_t kk=0;kk<nbOfSubEdges1 && !found;kk++)
 566                     {
 567                       found=(subEdge1PossiblyAlreadyIn1[2*kk]==idBg && subEdge1PossiblyAlreadyIn1[2*kk+1]==idEnd) || (subEdge1PossiblyAlreadyIn1[2*kk]==idEnd && subEdge1PossiblyAlreadyIn1[2*kk+1]==idBg);
 568                       if(!found)
 569                         offset2++;
 570                     }
 571                   if(found)
 572                     {
 573                       ElementaryEdge *e=pol1[offset1+(direct1?offset2:nbOfSubEdges1-offset2-1)];
 574                       e->getPtr()->declareOn();
 575                     }
 576                 }
 577             }
 578           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
 579         }
 580     }
 581 }
 582
 583 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
 584 {
 585   int nbOfNodesInPg=0;
 586   bool presenceOfQuadratic=presenceOfQuadraticEdge();
 587   conn.push_back(presenceOfQuadratic?NORM_QPOLYG:NORM_POLYGON);
 588   for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++)
 589     {
 590       Node *tmp=0;
 591       tmp=(*it)->getStartNode();
 592       std::map<INTERP_KERNEL::Node *,int>::const_iterator it1=mapp.find(tmp);
 593       conn.push_back((*it1).second);
 594       nbOfNodesInPg++;
 595     }
 596   if(presenceOfQuadratic)
 597     {
 598       int j=0;
 599       int off=offset+((int)addCoordsQuadratic.size())/2;
 600       for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++,j++,nbOfNodesInPg++)
 601         {
 602           INTERP_KERNEL::Node *node=(*it)->getPtr()->buildRepresentantOfMySelf();
 603           node->unApplySimilarity(xBary,yBary,fact);
 604           addCoordsQuadratic.push_back((*node)[0]);
 605           addCoordsQuadratic.push_back((*node)[1]);
 606           conn.push_back(off+j);
 607           node->decrRef();
 608         }
 609     }
 610   connI.push_back(connI.back()+nbOfNodesInPg+1);
 611 }
 612
 613 /*!
 614  * This method make the hypothesis that 'this' and 'other' are splited at the minimum into edges that are fully IN, OUT or ON.
 615  * This method returns newly created polygons in 'conn' and 'connI' and the corresponding ids ('idThis','idOther') are stored respectively into 'nbThis' and 'nbOther'.
 616  * @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
 617  * @param [in,out] edgesBoundaryOther, parameter that strores all edges in result of intersection that are not
 618  */
 619 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)
 620 {
 621   double xBaryBB, yBaryBB;
 622   double fact=normalizeExt(&other, xBaryBB, yBaryBB);
 623   //Locate 'this' relative to 'other'
 624   other.performLocatingOperationSlow(*this);  // without any assumption
 625   std::vector<QuadraticPolygon *> res=buildIntersectionPolygons(other,*this);
 626   for(std::vector<QuadraticPolygon *>::iterator it=res.begin();it!=res.end();it++)
 627     {
 628       (*it)->appendCrudeData(mapp,xBaryBB,yBaryBB,fact,offset,addCoordsQuadratic,conn,connI);
 629       INTERP_KERNEL::IteratorOnComposedEdge it1(*it);
 630       for(it1.first();!it1.finished();it1.next())
 631         {
 632           Edge *e=it1.current()->getPtr();
 633           if(edgesThis.find(e)!=edgesThis.end())
 634             edgesThis.erase(e);
 635           else
 636             {
 637               if(edgesBoundaryOther.find(e)!=edgesBoundaryOther.end())
 638                 edgesBoundaryOther.erase(e);
 639               else
 640                 edgesBoundaryOther.insert(e);
 641             }
 642         }
 643       nbThis.push_back(idThis);
 644       nbOther.push_back(idOther);
 645       delete *it;
 646     }
 647   unApplyGlobalSimilarityExt(other,xBaryBB,yBaryBB,fact);
 648 }
 649
 650 /*!
 651  * Warning This method is \b NOT const. 'this' and 'other' are modified after call of this method.
 652  * 'other' is a QuadraticPolygon of \b non closed edges.
 653  */
 654 double QuadraticPolygon::intersectWithAbs1D(QuadraticPolygon& other, bool& isColinear)
 655 {
 656   double ret = 0., xBaryBB, yBaryBB;
 657   double fact = normalize(&other, xBaryBB, yBaryBB);
 658
 659   QuadraticPolygon cpyOfThis(*this);
 660   QuadraticPolygon cpyOfOther(other);
 661   int nbOfSplits = 0;
 662   SplitPolygonsEachOther(cpyOfThis, cpyOfOther, nbOfSplits);
 663   //At this point cpyOfThis and cpyOfOther have been splited at maximum edge so that in/out can been done.
 664   performLocatingOperation(cpyOfOther);
 665   isColinear = false;
 666   for(std::list<ElementaryEdge *>::const_iterator it=cpyOfOther._sub_edges.begin();it!=cpyOfOther._sub_edges.end();it++)
 667     {
 668       switch((*it)->getLoc())
 669         {
 670         case FULL_IN_1:
 671           {
 672             ret += fabs((*it)->getPtr()->getCurveLength());
 673             break;
 674           }
 675         case FULL_ON_1:
 676           {
 677             isColinear=true;
 678             ret += fabs((*it)->getPtr()->getCurveLength());
 679             break;
 680           }
 681         default:
 682           {
 683           }
 684         }
 685     }
 686   return ret * fact;
 687 }
 688
 689 /*!
 690  * Warning contrary to intersectWith method this method is \b NOT const. 'this' and 'other' are modified after call of this method.
 691  */
 692 double QuadraticPolygon::intersectWithAbs(QuadraticPolygon& other, double* barycenter)
 693 {
 694   double ret=0.,bary[2],area,xBaryBB,yBaryBB;
 695   barycenter[0] = barycenter[1] = 0.;
 696   double fact=normalize(&other,xBaryBB,yBaryBB);
 697   std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
 698   for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
 699     {
 700       area=fabs((*iter)->getArea());
 701       (*iter)->getBarycenter(bary);
 702       delete *iter;
 703       ret+=area;
 704       barycenter[0] += bary[0]*area;
 705       barycenter[1] += bary[1]*area;
 706     }
 707   if ( ret > std::numeric_limits<double>::min() )
 708     {
 709       barycenter[0]=barycenter[0]/ret*fact+xBaryBB;
 710       barycenter[1]=barycenter[1]/ret*fact+yBaryBB;
 711
 712     }
 713   return ret*fact*fact;
 714 }
 715
 716 /*!
 717  * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
 718  * This is possible because loc attribute in Edge class is mutable.
 719  * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
 720  */
 721 double QuadraticPolygon::intersectWith(const QuadraticPolygon& other) const
 722 {
 723   double ret=0.;
 724   std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
 725   for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
 726     {
 727       ret+=fabs((*iter)->getArea());
 728       delete *iter;
 729     }
 730   return ret;
 731 }
 732
 733 /*!
 734  * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
 735  * This is possible because loc attribute in Edge class is mutable.
 736  * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
 737  */
 738 double QuadraticPolygon::intersectWith(const QuadraticPolygon& other, double* barycenter) const
 739 {
 740   double ret=0., bary[2];
 741   barycenter[0] = barycenter[1] = 0.;
 742   std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
 743   for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
 744     {
 745       double area = fabs((*iter)->getArea());
 746       (*iter)->getBarycenter(bary);
 747       delete *iter;
 748       ret+=area;
 749       barycenter[0] += bary[0]*area;
 750       barycenter[1] += bary[1]*area;
 751     }
 752   if ( ret > std::numeric_limits<double>::min() )
 753     {
 754       barycenter[0] /= ret;
 755       barycenter[1] /= ret;
 756     }
 757   return ret;
 758 }
 759
 760 /*!
 761  * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
 762  * This is possible because loc attribute in Edge class is mutable.
 763  * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
 764  */
 765 void QuadraticPolygon::intersectForPerimeter(const QuadraticPolygon& other, double& perimeterThisPart, double& perimeterOtherPart, double& perimeterCommonPart) const
 766 {
 767   perimeterThisPart=0.; perimeterOtherPart=0.; perimeterCommonPart=0.;
 768   QuadraticPolygon cpyOfThis(*this);
 769   QuadraticPolygon cpyOfOther(other); int nbOfSplits=0;
 770   SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits);
 771   performLocatingOperation(cpyOfOther);
 772   other.performLocatingOperation(cpyOfThis);
 773   cpyOfThis.dispatchPerimeterExcl(perimeterThisPart,perimeterCommonPart);
 774   cpyOfOther.dispatchPerimeterExcl(perimeterOtherPart,perimeterCommonPart);
 775   perimeterCommonPart/=2.;
 776 }
 777
 778 /*!
 779  * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
 780  * This is possible because loc attribute in Edge class is mutable.
 781  * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
 782  *
 783  * polThis.size()==this->size() and polOther.size()==other.size().
 784  * For each ElementaryEdge of 'this', the corresponding contribution in resulting polygon is in 'polThis'.
 785  * For each ElementaryEdge of 'other', the corresponding contribution in resulting polygon is in 'polOther'.
 786  * As consequence common part are counted twice (in polThis \b and in polOther).
 787  */
 788 void QuadraticPolygon::intersectForPerimeterAdvanced(const QuadraticPolygon& other, std::vector< double >& polThis, std::vector< double >& polOther) const
 789 {
 790   polThis.resize(size());
 791   polOther.resize(other.size());
 792   IteratorOnComposedEdge it1(const_cast<QuadraticPolygon *>(this));
 793   int edgeId=0;
 794   for(it1.first();!it1.finished();it1.next(),edgeId++)
 795     {
 796       ElementaryEdge* curE1=it1.current();
 797       QuadraticPolygon cpyOfOther(other);
 798       QuadraticPolygon tmp;
 799       tmp.pushBack(curE1->clone());
 800       int tmp2;
 801       SplitPolygonsEachOther(tmp,cpyOfOther,tmp2);
 802       other.performLocatingOperation(tmp);
 803       tmp.dispatchPerimeter(polThis[edgeId]);
 804     }
 805   //
 806   IteratorOnComposedEdge it2(const_cast<QuadraticPolygon *>(&other));
 807   edgeId=0;
 808   for(it2.first();!it2.finished();it2.next(),edgeId++)
 809     {
 810       ElementaryEdge* curE2=it2.current();
 811       QuadraticPolygon cpyOfThis(*this);
 812       QuadraticPolygon tmp;
 813       tmp.pushBack(curE2->clone());
 814       int tmp2;
 815       SplitPolygonsEachOther(tmp,cpyOfThis,tmp2);
 816       performLocatingOperation(tmp);
 817       tmp.dispatchPerimeter(polOther[edgeId]);
 818     }
 819 }
 820
 821
 822 /*!
 823  * numberOfCreatedPointsPerEdge is resized to the number of edges of 'this'.
 824  * This method returns in ordered maner the number of newly created points per edge.
 825  * This method performs a split process between 'this' and 'other' that gives the result PThis.
 826  * Then for each edges of 'this' this method counts how many edges in Pthis have the same id.
 827  */
 828 void QuadraticPolygon::intersectForPoint(const QuadraticPolygon& other, std::vector< int >& numberOfCreatedPointsPerEdge) const
 829 {
 830   numberOfCreatedPointsPerEdge.resize(size());
 831   IteratorOnComposedEdge it1(const_cast<QuadraticPolygon *>(this));
 832   int edgeId=0;
 833   for(it1.first();!it1.finished();it1.next(),edgeId++)
 834     {
 835       ElementaryEdge* curE1=it1.current();
 836       QuadraticPolygon cpyOfOther(other);
 837       QuadraticPolygon tmp;
 838       tmp.pushBack(curE1->clone());
 839       int tmp2;
 840       SplitPolygonsEachOther(tmp,cpyOfOther,tmp2);
 841       numberOfCreatedPointsPerEdge[edgeId]=tmp.recursiveSize()-1;
 842     }
 843 }
 844
 845 /*!
 846  * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
 847  * This is possible because loc attribute in Edge class is mutable.
 848  * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
 849  */
 850 std::vector<QuadraticPolygon *> QuadraticPolygon::intersectMySelfWith(const QuadraticPolygon& other) const
 851 {
 852   QuadraticPolygon cpyOfThis(*this);
 853   QuadraticPolygon cpyOfOther(other); int nbOfSplits=0;
 854   SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits);
 855   //At this point cpyOfThis and cpyOfOther have been splited at maximum edge so that in/out can been done.
 856   performLocatingOperation(cpyOfOther);
 857   return other.buildIntersectionPolygons(cpyOfThis,cpyOfOther);
 858 }
 859
 860 /*!
 861  * This method is typically the first step of boolean operations between pol1 and pol2.
 862  * This method perform the minimal splitting so that at the end each edges constituting pol1 are fully either IN or OUT or ON.
 863  * @param pol1 IN/OUT param that is equal to 'this' when called.
 864  */
 865 void QuadraticPolygon::SplitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticPolygon& pol2, int& nbOfSplits)
 866 {
 867   IteratorOnComposedEdge it1(&pol1),it2(&pol2);
 868   MergePoints merge;
 869   ComposedEdge *c1=new ComposedEdge;
 870   ComposedEdge *c2=new ComposedEdge;
 871   for(it2.first();!it2.finished();it2.next())
 872     {
 873       ElementaryEdge* curE2=it2.current();
 874       if(!curE2->isThereStartPoint())
 875         it1.first();
 876       else
 877         it1=curE2->getIterator();
 878       for(;!it1.finished();)
 879         {
 880
 881           ElementaryEdge* curE1=it1.current();
 882           merge.clear(); nbOfSplits++;
 883           if(curE1->getPtr()->intersectWith(curE2->getPtr(),merge,*c1,*c2))
 884             {
 885               if(!curE1->getDirection()) c1->reverse();
 886               if(!curE2->getDirection()) c2->reverse();
 887               UpdateNeighbours(merge,it1,it2,c1,c2);
 888               //Substitution of simple edge by sub-edges.
 889               delete curE1; // <-- destroying simple edge coming from pol1
 890               delete curE2; // <-- destroying simple edge coming from pol2
 891               it1.insertElemEdges(c1,true);// <-- 2nd param is true to go next.
 892               it2.insertElemEdges(c2,false);// <-- 2nd param is false to avoid to go next.
 893               curE2=it2.current();
 894               //
 895               it1.assignMySelfToAllElems(c2);//To avoid that others
 896               SoftDelete(c1);
 897               SoftDelete(c2);
 898               c1=new ComposedEdge;
 899               c2=new ComposedEdge;
 900             }
 901           else
 902             {
 903               UpdateNeighbours(merge,it1,it2,curE1,curE2);
 904               it1.next();
 905             }
 906         }
 907     }
 908   Delete(c1);
 909   Delete(c2);
 910 }
 911
 912 void QuadraticPolygon::performLocatingOperation(QuadraticPolygon& pol2) const
 913 {
 914   IteratorOnComposedEdge it(&pol2);
 915   TypeOfEdgeLocInPolygon loc=FULL_ON_1;
 916   for(it.first();!it.finished();it.next())
 917     {
 918       ElementaryEdge *cur=it.current();
 919       loc=cur->locateFullyMySelf(*this,loc);
 920     }
 921 }
 922
 923 void QuadraticPolygon::performLocatingOperationSlow(QuadraticPolygon& pol2) const
 924 {
 925   IteratorOnComposedEdge it(&pol2);
 926   for(it.first();!it.finished();it.next())
 927     {
 928       ElementaryEdge *cur=it.current();
 929       cur->locateFullyMySelfAbsolute(*this);
 930     }
 931 }
 932
 933 /*!
 934  * Given 2 polygons 'pol1' and 'pol2' (localized) the resulting polygons are returned.
 935  *
 936  * this : pol2 simplified.
 937  * @param pol1 pol1 split.
 938  * @param pol2 pol2 split.
 939  */
 940 std::vector<QuadraticPolygon *> QuadraticPolygon::buildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2) const
 941 {
 942   std::vector<QuadraticPolygon *> ret;
 943   std::list<QuadraticPolygon *> pol2Zip=pol2.zipConsecutiveInSegments();
 944   if(!pol2Zip.empty())
 945     closePolygons(pol2Zip,pol1,ret);
 946   else
 947     {//borders of pol2 do not cross pol1,and pol2 borders are outside of pol1. That is to say, either pol2 and pol1
 948       //do not overlap or  pol1 is fully inside pol2. So in the first case no intersection, in the other case
 949       //the intersection is pol1.
 950       ElementaryEdge *e1FromPol1=pol1[0];
 951       TypeOfEdgeLocInPolygon loc=FULL_ON_1;
 952       loc=e1FromPol1->locateFullyMySelf(*this,loc);
 953       if(loc==FULL_IN_1)
 954         ret.push_back(new QuadraticPolygon(pol1));
 955     }
 956   return ret;
 957 }
 958
 959 /*!
 960  * Returns parts of potentially non closed-polygons. Each returned polygons are not mergeable.
 961  * this : pol2 split and locallized.
 962  */
 963 std::list<QuadraticPolygon *> QuadraticPolygon::zipConsecutiveInSegments() const
 964 {
 965   std::list<QuadraticPolygon *> ret;
 966   IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(this));
 967   int nbOfTurns=recursiveSize();
 968   int i=0;
 969   if(!it.goToNextInOn(false,i,nbOfTurns))
 970     return ret;
 971   i=0;
 972   //
 973   while(i<nbOfTurns)
 974     {
 975       QuadraticPolygon *tmp1=new QuadraticPolygon;
 976       TypeOfEdgeLocInPolygon loc=it.current()->getLoc();
 977       while(loc!=FULL_OUT_1 && i<nbOfTurns)
 978         {
 979           ElementaryEdge *tmp3=it.current()->clone();
 980           tmp1->pushBack(tmp3);
 981           it.nextLoop(); i++;
 982           loc=it.current()->getLoc();
 983         }
 984       if(tmp1->empty())
 985         {
 986           delete tmp1;
 987           continue;
 988         }
 989       ret.push_back(tmp1);
 990       it.goToNextInOn(true,i,nbOfTurns);
 991     }
 992   return ret;
 993 }
 994
 995 /*!
 996  * 'this' should be considered as pol2Simplified.
 997  * @param pol2zip is a list of set of edges (openned polygon) coming from split polygon 2.
 998  * @param pol1 is split pol1.
 999  * @param results the resulting \b CLOSED polygons.
1000  */
1001 void QuadraticPolygon::closePolygons(std::list<QuadraticPolygon *>& pol2Zip, const QuadraticPolygon& pol1,
1002                                      std::vector<QuadraticPolygon *>& results) const
1003 {
1004   bool directionKnownInPol1=false;
1005   bool directionInPol1;
1006   for(std::list<QuadraticPolygon *>::iterator iter=pol2Zip.begin();iter!=pol2Zip.end();)
1007     {
1008       if((*iter)->completed())
1009         {
1010           results.push_back(*iter);
1011           directionKnownInPol1=false;
1012           iter=pol2Zip.erase(iter);
1013           continue;
1014         }
1015       if(!directionKnownInPol1)
1016         {
1017           if(!(*iter)->amIAChanceToBeCompletedBy(pol1,*this,directionInPol1))
1018             { delete *iter; iter=pol2Zip.erase(iter); continue; }
1019           else
1020             directionKnownInPol1=true;
1021         }
1022       std::list<QuadraticPolygon *>::iterator iter2=iter; iter2++;
1023       std::list<QuadraticPolygon *>::iterator iter3=(*iter)->fillAsMuchAsPossibleWith(pol1,iter2,pol2Zip.end(),directionInPol1);
1024       if(iter3!=pol2Zip.end())
1025         {
1026           (*iter)->pushBack(*iter3);
1027           SoftDelete(*iter3);
1028           pol2Zip.erase(iter3);
1029         }
1030     }
1031 }
1032
1033 /*!
1034  * 'this' is expected to be set of edges (not closed) of pol2 split.
1035  */
1036 bool QuadraticPolygon::amIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted,const QuadraticPolygon& pol2NotSplitted, bool& direction)
1037 {
1038   IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(&pol1Splitted));
1039   bool found=false;
1040   Node *n=getEndNode();
1041   ElementaryEdge *cur=it.current();
1042   for(it.first();!it.finished() && !found;)
1043     {
1044       cur=it.current();
1045       found=(cur->getStartNode()==n);
1046       if(!found)
1047         it.next();
1048     }
1049   if(!found)
1050     throw Exception("Internal error : polygons uncompatible each others. Should never happend");
1051   //Ok we found correspondance between this and pol1. Searching for right direction to close polygon.
1052   ElementaryEdge *e=_sub_edges.back();
1053   if(e->getLoc()==FULL_ON_1)
1054     {
1055       if(e->getPtr()==cur->getPtr())
1056         {
1057           direction=false;
1058           it.previousLoop();
1059           cur=it.current();
1060           Node *repr=cur->getPtr()->buildRepresentantOfMySelf();
1061           bool ret=pol2NotSplitted.isInOrOut(repr);
1062           repr->decrRef();
1063           return ret;
1064         }
1065       else
1066         {
1067           direction=true;
1068           Node *repr=cur->getPtr()->buildRepresentantOfMySelf();
1069           bool ret=pol2NotSplitted.isInOrOut(repr);
1070           repr->decrRef();
1071           return ret;
1072         }
1073     }
1074   else
1075     direction=cur->locateFullyMySelfAbsolute(pol2NotSplitted)==FULL_IN_1;
1076   return true;
1077 }
1078
1079 /*!
1080  * This method fills as much as possible 'this' (part of pol2 split) with edges of 'pol1Splitted'.
1081  */
1082 std::list<QuadraticPolygon *>::iterator QuadraticPolygon::fillAsMuchAsPossibleWith(const QuadraticPolygon& pol1Splitted,
1083                                                                                    std::list<QuadraticPolygon *>::iterator iStart,
1084                                                                                    std::list<QuadraticPolygon *>::iterator iEnd,
1085                                                                                    bool direction)
1086 {
1087   IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(&pol1Splitted));
1088   bool found=false;
1089   Node *n=getEndNode();
1090   ElementaryEdge *cur;
1091   for(it.first();!it.finished() && !found;)
1092     {
1093       cur=it.current();
1094       found=(cur->getStartNode()==n);
1095       if(!found)
1096         it.next();
1097     }
1098   if(!direction)
1099     it.previousLoop();
1100   Node *nodeToTest;
1101   std::list<QuadraticPolygon *>::iterator ret;
1102   do
1103     {
1104       cur=it.current();
1105       ElementaryEdge *tmp=cur->clone();
1106       if(!direction)
1107         tmp->reverse();
1108       pushBack(tmp);
1109       nodeToTest=tmp->getEndNode();
1110       direction?it.nextLoop():it.previousLoop();
1111       ret=CheckInList(nodeToTest,iStart,iEnd);
1112       if(completed())
1113         return iEnd;
1114     }
1115   while(ret==iEnd);
1116   return ret;
1117 }
1118
1119 std::list<QuadraticPolygon *>::iterator QuadraticPolygon::CheckInList(Node *n, std::list<QuadraticPolygon *>::iterator iStart,
1120                                                                       std::list<QuadraticPolygon *>::iterator iEnd)
1121 {
1122   for(std::list<QuadraticPolygon *>::iterator iter=iStart;iter!=iEnd;iter++)
1123     if((*iter)->isNodeIn(n))
1124       return iter;
1125   return iEnd;
1126 }
1127
1128 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,
1129                                        std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI, std::vector<int>& nb1, std::vector<int>& nb2)
1130 {
1131   pol1.initLocations();
1132   for(std::set<Edge *>::const_iterator it9=notUsedInPol1.begin();it9!=notUsedInPol1.end();it9++)
1133     { (*it9)->initLocs(); (*it9)->declareOn(); }
1134   for(std::set<Edge *>::const_iterator itA=edgesInPol2OnBoundary.begin();itA!=edgesInPol2OnBoundary.end();itA++)
1135     { (*itA)->initLocs(); (*itA)->declareIn(); }
1136   ////
1137   std::set<Edge *> notUsedInPol1L(notUsedInPol1);
1138   IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(&pol1));
1139   int sz=pol1.size();
1140   std::list<QuadraticPolygon *> pol1Zip;
1141   if(pol1.size()==(int)notUsedInPol1.size() && edgesInPol2OnBoundary.empty())
1142     {
1143       pol1.appendCrudeData(mapp,0.,0.,1.,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1);
1144       return ;
1145     }
1146   while(!notUsedInPol1L.empty())
1147     {
1148       for(int i=0;i<sz && (it.current()->getStartNode()->getLoc()!=IN_1 || it.current()->getLoc()!=FULL_ON_1);i++)
1149         it.nextLoop();
1150       if(it.current()->getStartNode()->getLoc()!=IN_1 || it.current()->getLoc()!=FULL_ON_1)
1151         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 !");
1152       QuadraticPolygon *tmp1=new QuadraticPolygon;
1153       do
1154         {
1155           Edge *ee=it.current()->getPtr();
1156           if(ee->getLoc()==FULL_ON_1)
1157             {
1158               ee->incrRef(); notUsedInPol1L.erase(ee);
1159               tmp1->pushBack(new ElementaryEdge(ee,it.current()->getDirection()));
1160             }
1161           it.nextLoop();
1162         }
1163       while(it.current()->getStartNode()->getLoc()!=IN_1 && !notUsedInPol1L.empty());
1164       pol1Zip.push_back(tmp1);
1165     }
1166   ////
1167   std::list<QuadraticPolygon *> retPolsUnderContruction;
1168   std::list<Edge *> edgesInPol2OnBoundaryL(edgesInPol2OnBoundary.begin(),edgesInPol2OnBoundary.end());
1169   std::map<QuadraticPolygon *, std::list<QuadraticPolygon *> > pol1ZipConsumed;
1170   std::size_t maxNbOfTurn=edgesInPol2OnBoundaryL.size(),nbOfTurn=0,iiMNT=0;
1171   for(std::list<QuadraticPolygon *>::const_iterator itMNT=pol1Zip.begin();itMNT!=pol1Zip.end();itMNT++,iiMNT++)
1172     nbOfTurn+=(*itMNT)->size();
1173   maxNbOfTurn=maxNbOfTurn*nbOfTurn; maxNbOfTurn*=maxNbOfTurn;
1174   nbOfTurn=0;
1175   while(nbOfTurn<maxNbOfTurn && ((!pol1Zip.empty() || !edgesInPol2OnBoundaryL.empty())))
1176     {
1177       for(std::list<QuadraticPolygon *>::iterator it1=retPolsUnderContruction.begin();it1!=retPolsUnderContruction.end();)
1178         {
1179           if((*it1)->getStartNode()==(*it1)->getEndNode())
1180             {
1181               it1++;
1182               continue;
1183             }
1184           Node *curN=(*it1)->getEndNode();
1185           bool smthHappened=false;
1186           for(std::list<Edge *>::iterator it2=edgesInPol2OnBoundaryL.begin();it2!=edgesInPol2OnBoundaryL.end();)
1187             {
1188               if(curN==(*it2)->getStartNode())
1189                 { (*it2)->incrRef(); (*it1)->pushBack(new ElementaryEdge(*it2,true)); curN=(*it2)->getEndNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); }
1190               else if(curN==(*it2)->getEndNode())
1191                 { (*it2)->incrRef(); (*it1)->pushBack(new ElementaryEdge(*it2,false)); curN=(*it2)->getStartNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); }
1192               else
1193                 it2++;
1194             }
1195           if(smthHappened)
1196             {
1197               for(std::list<QuadraticPolygon *>::iterator it3=pol1Zip.begin();it3!=pol1Zip.end();)
1198                 {
1199                   if(curN==(*it3)->getStartNode())
1200                     {
1201                       for(std::list<ElementaryEdge *>::const_iterator it4=(*it3)->_sub_edges.begin();it4!=(*it3)->_sub_edges.end();it4++)
1202                         { (*it4)->getPtr()->incrRef(); bool dir=(*it4)->getDirection(); (*it1)->pushBack(new ElementaryEdge((*it4)->getPtr(),dir)); }
1203                       smthHappened=true;
1204                       pol1ZipConsumed[*it1].push_back(*it3);
1205                       curN=(*it3)->getEndNode();
1206                       it3=pol1Zip.erase(it3);
1207                     }
1208                   else
1209                     it3++;
1210                 }
1211             }
1212           if(!smthHappened)
1213             {
1214               for(std::list<ElementaryEdge *>::const_iterator it5=(*it1)->_sub_edges.begin();it5!=(*it1)->_sub_edges.end();it5++)
1215                 {
1216                   Edge *ee=(*it5)->getPtr();
1217                   if(edgesInPol2OnBoundary.find(ee)!=edgesInPol2OnBoundary.end())
1218                     edgesInPol2OnBoundaryL.push_back(ee);
1219                 }
1220               for(std::list<QuadraticPolygon *>::iterator it6=pol1ZipConsumed[*it1].begin();it6!=pol1ZipConsumed[*it1].end();it6++)
1221                 pol1Zip.push_front(*it6);
1222               pol1ZipConsumed.erase(*it1);
1223               delete *it1;
1224               it1=retPolsUnderContruction.erase(it1);
1225             }
1226         }
1227       if(!pol1Zip.empty())
1228         {
1229           QuadraticPolygon *tmp=new QuadraticPolygon;
1230           QuadraticPolygon *first=*(pol1Zip.begin());
1231           for(std::list<ElementaryEdge *>::const_iterator it4=first->_sub_edges.begin();it4!=first->_sub_edges.end();it4++)
1232             { (*it4)->getPtr()->incrRef(); bool dir=(*it4)->getDirection(); tmp->pushBack(new ElementaryEdge((*it4)->getPtr(),dir)); }
1233           pol1ZipConsumed[tmp].push_back(first);
1234           retPolsUnderContruction.push_back(tmp);
1235           pol1Zip.erase(pol1Zip.begin());
1236         }
1237       nbOfTurn++;
1238     }
1239   if(nbOfTurn==maxNbOfTurn)
1240     {
1241       std::ostringstream oss; oss << "Error during reconstruction of residual of cell ! It appears that either source or/and target mesh is/are not conform !";
1242       oss << " Number of turns is = " << nbOfTurn << " !";
1243       throw INTERP_KERNEL::Exception(oss.str().c_str());
1244     }
1245   for(std::list<QuadraticPolygon *>::iterator it1=retPolsUnderContruction.begin();it1!=retPolsUnderContruction.end();it1++)
1246     {
1247       if((*it1)->getStartNode()==(*it1)->getEndNode())
1248         {
1249           (*it1)->appendCrudeData(mapp,0.,0.,1.,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1);
1250           for(std::list<QuadraticPolygon *>::iterator it6=pol1ZipConsumed[*it1].begin();it6!=pol1ZipConsumed[*it1].end();it6++)
1251             delete *it6;
1252           delete *it1;
1253         }
1254     }
1255 }