src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx

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