src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx

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