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