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 itThis(this),itOther(&other);  // other is (part of) the tool mesh
 299   MergePoints merge;
 300   ComposedEdge *cThis=new ComposedEdge;
 301   ComposedEdge *cOther=new ComposedEdge;
 302   int i=0;
 303   std::map<INTERP_KERNEL::Node *,int> mapAddCoo;
 304   for(itOther.first();!itOther.finished();itOther.next(),i++)
 305     {
 306       // For each edge of 'other', proceed with intersections: the edge might split into sub-edges, 'otherTmp' will hold the final split result.
 307       // In the process of going through all edges of 'other', 'this' (which contains initially only one edge)
 308       // is sub-divised into several edges : each of them has to be tested when intersecting the next candidate stored in 'other'.
 309       QuadraticPolygon otherTmp;
 310       ElementaryEdge* curOther=itOther.current();
 311       otherTmp.pushBack(new ElementaryEdge(curOther->getPtr(),curOther->getDirection())); curOther->getPtr()->incrRef();
 312       IteratorOnComposedEdge itOtherTmp(&otherTmp);
 313       for(itOtherTmp.first();!itOtherTmp.finished();itOtherTmp.next())
 314         {
 315           ElementaryEdge* curOtherTmp=itOtherTmp.current();
 316           if(!curOtherTmp->isThereStartPoint())
 317             itThis.first();   // reset iterator on 'this'
 318           else
 319             itThis=curOtherTmp->getIterator();
 320           for(;!itThis.finished();)
 321             {
 322               ElementaryEdge* curThis=itThis.current();
 323               merge.clear();
 324               //
 325               std::map<INTERP_KERNEL::Node *,int>::const_iterator thisStart(mapThis.find(curThis->getStartNode())),thisEnd(mapThis.find(curThis->getEndNode())),
 326                                                                   otherStart(mapOther.find(curOtherTmp->getStartNode())),otherEnd(mapOther.find(curOtherTmp->getEndNode()));
 327               int thisStart2(thisStart==mapThis.end()?-1:(*thisStart).second), thisEnd2(thisEnd==mapThis.end()?-1:(*thisEnd).second),
 328                   otherStart2(otherStart==mapOther.end()?-1:(*otherStart).second+offset1),otherEnd2(otherEnd==mapOther.end()?-1:(*otherEnd).second+offset1);
 329               //
 330               if(curThis->getPtr()->intersectWith(curOtherTmp->getPtr(),merge,*cThis,*cOther))
 331                 {
 332                   if(!curThis->getDirection()) cThis->reverse();
 333                   if(!curOtherTmp->getDirection()) cOther->reverse();
 334                   // Substitution of a single simple edge by two sub-edges resulting from the intersection
 335                   //    First modify the edges currently pointed by itThis and itOtherTmp so that the newly created node
 336                   //    becomes the end of the previous sub-edge and the beginning of the next one.
 337                   UpdateNeighbours(merge,itThis,itOtherTmp,cThis,cOther);
 338                   delete curThis;           // <-- destroying simple edge coming from pol1
 339                   delete curOtherTmp;       // <-- destroying simple edge coming from pol2
 340                   //     Then insert second part of the intersection.
 341                   itThis.insertElemEdges(cThis,true);       // <-- 2nd param is true to go next.
 342                   itOtherTmp.insertElemEdges(cOther,false); // <-- 2nd param is false to avoid to go next.
 343                   curOtherTmp=itOtherTmp.current();
 344                   //
 345                   itThis.assignMySelfToAllElems(cOther);
 346                   SoftDelete(cThis);
 347                   SoftDelete(cOther);
 348                   cThis=new ComposedEdge;
 349                   cOther=new ComposedEdge;
 350                 }
 351               else
 352                 {
 353                   UpdateNeighbours(merge,itThis,itOtherTmp,curThis,curOtherTmp);
 354                   itThis.next();
 355                 }
 356               merge.updateMergedNodes(thisStart2,thisEnd2,otherStart2,otherEnd2,mergedNodes);
 357             }
 358         }
 359       // If one sub-edge of otherTmp is "ON" an edge of this, then we have colinearity (all edges in otherTmp are //)
 360       if(otherTmp.presenceOfOn())
 361         edgesInOtherColinearWithThis[otherEdgeIds[i]].push_back(cellIdThis);
 362       // Converting back to integer connectivity:
 363       if(otherTmp._sub_edges.size()>1)   // only if a new point has been added (i.e. an actual intersection was done)
 364         {
 365           for(std::list<ElementaryEdge *>::const_iterator it=otherTmp._sub_edges.begin();it!=otherTmp._sub_edges.end();it++)
 366             (*it)->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,/**/fact,xBaryBB,yBaryBB,/**/subDivOther[otherEdgeIds[i]],addCoo,mapAddCoo);
 367         }
 368     }
 369   Delete(cThis);
 370   Delete(cOther);
 371   //
 372   for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++)
 373     (*it)->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,/**/fact,xBaryBB,yBaryBB,/**/edgesThis,addCoo,mapAddCoo);
 374   //
 375 }
 376
 377 /*!
 378  * This method builds 'this' from its descending conn stored in crude mode (MEDCoupling).
 379  * Descending conn is in FORTRAN relative mode in order to give the
 380  * orientation of edge (see buildDescendingConnectivity2() method).
 381  * See appendEdgeFromCrudeDataArray() for params description.
 382  */
 383 void QuadraticPolygon::buildFromCrudeDataArray(const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad, const int *nodalBg, const double *coords,
 384                                                const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges)
 385 {
 386   std::size_t nbOfSeg=std::distance(descBg,descEnd);
 387   for(std::size_t i=0;i<nbOfSeg;i++)
 388     {
 389       appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
 390     }
 391 }
 392
 393 void QuadraticPolygon::appendEdgeFromCrudeDataArray(std::size_t edgePos, const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad,
 394                                                     const int *nodalBg, const double *coords,
 395                                                     const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges)
 396 {
 397   if(!isQuad)
 398     {
 399       bool direct=descBg[edgePos]>0;
 400       int edgeId=abs(descBg[edgePos])-1; // back to C indexing mode
 401       const std::vector<int>& subEdge=intersectEdges[edgeId];
 402       std::size_t nbOfSubEdges=subEdge.size()/2;
 403       for(std::size_t j=0;j<nbOfSubEdges;j++)
 404         appendSubEdgeFromCrudeDataArray(0,j,direct,edgeId,subEdge,mapp);
 405     }
 406   else
 407     {
 408       std::size_t nbOfSeg=std::distance(descBg,descEnd);
 409       const double *st=coords+2*(nodalBg[edgePos]);
 410       INTERP_KERNEL::Node *st0=new INTERP_KERNEL::Node(st[0],st[1]);
 411       const double *endd=coords+2*(nodalBg[(edgePos+1)%nbOfSeg]);
 412       INTERP_KERNEL::Node *endd0=new INTERP_KERNEL::Node(endd[0],endd[1]);
 413       const double *middle=coords+2*(nodalBg[edgePos+nbOfSeg]);
 414       INTERP_KERNEL::Node *middle0=new INTERP_KERNEL::Node(middle[0],middle[1]);
 415       EdgeLin *e1,*e2;
 416       e1=new EdgeLin(st0,middle0);
 417       e2=new EdgeLin(middle0,endd0);
 418       SegSegIntersector inters(*e1,*e2);
 419       bool colinearity=inters.areColinears();
 420       delete e1; delete e2;
 421       //
 422       bool direct=descBg[edgePos]>0;
 423       int edgeId=abs(descBg[edgePos])-1;
 424       const std::vector<int>& subEdge=intersectEdges[edgeId];
 425       std::size_t nbOfSubEdges=subEdge.size()/2;
 426       if(colinearity)
 427         {
 428           for(std::size_t j=0;j<nbOfSubEdges;j++)
 429             appendSubEdgeFromCrudeDataArray(0,j,direct,edgeId,subEdge,mapp);
 430         }
 431       else
 432         {
 433           Edge *e=new EdgeArcCircle(st0,middle0,endd0,true);
 434           for(std::size_t j=0;j<nbOfSubEdges;j++)
 435             appendSubEdgeFromCrudeDataArray(e,j,direct,edgeId,subEdge,mapp);
 436           e->decrRef();
 437         }
 438       st0->decrRef(); endd0->decrRef(); middle0->decrRef();
 439     }
 440 }
 441
 442 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)
 443 {
 444   std::size_t nbOfSubEdges=subEdge.size()/2;
 445   if(!baseEdge)
 446     {//it is not a quadratic subedge
 447       Node *start=(*mapp.find(direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1])).second;
 448       Node *end=(*mapp.find(direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2])).second;
 449       ElementaryEdge *e=ElementaryEdge::BuildEdgeFromStartEndDir(true,start,end);
 450       pushBack(e);
 451     }
 452   else
 453     {//it is a quadratic subedge
 454       Node *start=(*mapp.find(direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1])).second;
 455       Node *end=(*mapp.find(direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2])).second;
 456       Edge *ee=baseEdge->buildEdgeLyingOnMe(start,end);
 457       ElementaryEdge *eee=new ElementaryEdge(ee,true);
 458       pushBack(eee);
 459     }
 460 }
 461
 462 /*!
 463  * 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
 464  * orientation of edge.
 465  */
 466 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> >& intersectEdges2,
 467                                                 const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, const std::vector<std::vector<int> >& intersectEdges1,
 468                                                 const std::vector< std::vector<int> >& colinear1,
 469                                                 std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> >& alreadyExistingIn2)
 470 {
 471   std::size_t nbOfSeg=std::distance(descBg,descEnd);
 472   for(std::size_t i=0;i<nbOfSeg;i++)//loop over all edges of pol2
 473     {
 474       bool direct=descBg[i]>0;
 475       int edgeId=abs(descBg[i])-1;//current edge id of pol2
 476       std::map<int,std::vector<INTERP_KERNEL::ElementaryEdge *> >::const_iterator it1=alreadyExistingIn2.find(descBg[i]),it2=alreadyExistingIn2.find(-descBg[i]);
 477       if(it1!=alreadyExistingIn2.end() || it2!=alreadyExistingIn2.end())
 478         {
 479           bool sameDir=(it1!=alreadyExistingIn2.end());
 480           const std::vector<INTERP_KERNEL::ElementaryEdge *>& edgesAlreadyBuilt=sameDir?(*it1).second:(*it2).second;
 481           if(sameDir)
 482             {
 483               for(std::vector<INTERP_KERNEL::ElementaryEdge *>::const_iterator it3=edgesAlreadyBuilt.begin();it3!=edgesAlreadyBuilt.end();it3++)
 484                 {
 485                   Edge *ee=(*it3)->getPtr(); ee->incrRef();
 486                   pushBack(new ElementaryEdge(ee,(*it3)->getDirection()));
 487                 }
 488             }
 489           else
 490             {
 491               for(std::vector<INTERP_KERNEL::ElementaryEdge *>::const_reverse_iterator it4=edgesAlreadyBuilt.rbegin();it4!=edgesAlreadyBuilt.rend();it4++)
 492                 {
 493                   Edge *ee=(*it4)->getPtr(); ee->incrRef();
 494                   pushBack(new ElementaryEdge(ee,!(*it4)->getDirection()));
 495                 }
 496             }
 497           continue;
 498         }
 499       bool directos=colinear1[edgeId].empty();
 500       std::vector<std::pair<int,std::pair<bool,int> > > idIns1;
 501       int offset1=0;
 502       if(!directos)
 503         {// if the current edge of pol2 has one or more colinear edges part into pol1
 504           const std::vector<int>& c=colinear1[edgeId];
 505           std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1);
 506           for(std::size_t j=0;j<nbOfEdgesIn1;j++)
 507             {
 508               int edgeId1=abs(descBg1[j])-1;
 509               if(std::find(c.begin(),c.end(),edgeId1)!=c.end())
 510                 {
 511                   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
 512                   //std::pair<edgeId1); direct1=descBg1[j]>0;
 513                 }
 514               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
 515             }
 516           directos=idIns1.empty();
 517         }
 518       if(directos)
 519         {//no subpart of edge 'edgeId' of pol2 is in pol1 so let's operate the same thing that QuadraticPolygon::buildFromCrudeDataArray method
 520           std::size_t oldSz=_sub_edges.size();
 521           appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges2);
 522           std::size_t newSz=_sub_edges.size();
 523           std::size_t zeSz=newSz-oldSz;
 524           alreadyExistingIn2[descBg[i]].resize(zeSz);
 525           std::list<ElementaryEdge *>::const_reverse_iterator it5=_sub_edges.rbegin();
 526           for(std::size_t p=0;p<zeSz;p++,it5++)
 527             alreadyExistingIn2[descBg[i]][zeSz-p-1]=*it5;
 528         }
 529       else
 530         {//there is subpart of edge 'edgeId' of pol2 inside pol1
 531           const std::vector<int>& subEdge=intersectEdges2[edgeId];
 532           std::size_t nbOfSubEdges=subEdge.size()/2;
 533           for(std::size_t j=0;j<nbOfSubEdges;j++)
 534             {
 535               int idBg=direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1];
 536               int idEnd=direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2];
 537               bool direction11,found=false;
 538               bool direct1;//store if needed the direction in 1
 539               int offset2;
 540               std::size_t nbOfSubEdges1;
 541               for(std::vector<std::pair<int,std::pair<bool,int> > >::const_iterator it=idIns1.begin();it!=idIns1.end() && !found;it++)
 542                 {
 543                   int idIn1=(*it).first;//store if needed the cell id in 1
 544                   direct1=(*it).second.first;
 545                   offset1=(*it).second.second;
 546                   const std::vector<int>& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1];
 547                   nbOfSubEdges1=subEdge1PossiblyAlreadyIn1.size()/2;
 548                   offset2=0;
 549                   for(std::size_t k=0;k<nbOfSubEdges1 && !found;k++)
 550                     {//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
 551                       if(subEdge1PossiblyAlreadyIn1[2*k]==idBg && subEdge1PossiblyAlreadyIn1[2*k+1]==idEnd)
 552                         { direction11=true; found=true; }
 553                       else if(subEdge1PossiblyAlreadyIn1[2*k]==idEnd && subEdge1PossiblyAlreadyIn1[2*k+1]==idBg)
 554                         { direction11=false; found=true; }
 555                       else
 556                         offset2++;
 557                     }
 558                 }
 559               if(!found)
 560                 {//the current subedge of edge 'edgeId' of pol2 is not a part of the colinear edge 'idIn1' of pol1 -> build new Edge instance
 561                   //appendEdgeFromCrudeDataArray(j,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
 562                   Node *start=(*mapp.find(idBg)).second;
 563                   Node *end=(*mapp.find(idEnd)).second;
 564                   ElementaryEdge *e=ElementaryEdge::BuildEdgeFromStartEndDir(true,start,end);
 565                   pushBack(e);
 566                   alreadyExistingIn2[descBg[i]].push_back(e);
 567                 }
 568               else
 569                 {//the current subedge of edge 'edgeId' of pol2 is part of the colinear edge 'idIn1' of pol1 -> reuse Edge instance of pol1
 570                   ElementaryEdge *e=pol1[offset1+(direct1?offset2:nbOfSubEdges1-offset2-1)];
 571                   Edge *ee=e->getPtr();
 572                   ee->incrRef();
 573                   ElementaryEdge *e2=new ElementaryEdge(ee,!(direct1^direction11));
 574                   pushBack(e2);
 575                   alreadyExistingIn2[descBg[i]].push_back(e2);
 576                 }
 577             }
 578         }
 579     }
 580 }
 581
 582 /*!
 583  * Method expected to be called on pol2. Every params not suffixed by numbered are supposed to refer to pol2 (this).
 584  * Method to find edges that are ON.
 585  */
 586 void QuadraticPolygon::updateLocOfEdgeFromCrudeDataArray2(const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges,
 587                                                           const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1,
 588                                                           const std::vector<std::vector<int> >& intersectEdges1, const std::vector< std::vector<int> >& colinear1) const
 589 {
 590   std::size_t nbOfSeg=std::distance(descBg,descEnd);
 591   for(std::size_t i=0;i<nbOfSeg;i++)//loop over all edges of pol2
 592     {
 593       bool direct=descBg[i]>0;
 594       int edgeId=abs(descBg[i])-1;//current edge id of pol2
 595       const std::vector<int>& c=colinear1[edgeId];
 596       if(c.empty())
 597         continue;
 598       const std::vector<int>& subEdge=intersectEdges[edgeId];
 599       std::size_t nbOfSubEdges=subEdge.size()/2;
 600       //
 601       std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1);
 602       int offset1=0;
 603       for(std::size_t j=0;j<nbOfEdgesIn1;j++)
 604         {
 605           int edgeId1=abs(descBg1[j])-1;
 606           if(std::find(c.begin(),c.end(),edgeId1)!=c.end())
 607             {
 608               for(std::size_t k=0;k<nbOfSubEdges;k++)
 609                 {
 610                   int idBg=direct?subEdge[2*k]:subEdge[2*nbOfSubEdges-2*k-1];
 611                   int idEnd=direct?subEdge[2*k+1]:subEdge[2*nbOfSubEdges-2*k-2];
 612                   int idIn1=edgeId1;
 613                   bool direct1=descBg1[j]>0;
 614                   const std::vector<int>& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1];
 615                   std::size_t nbOfSubEdges1=subEdge1PossiblyAlreadyIn1.size()/2;
 616                   int offset2=0;
 617                   bool found=false;
 618                   for(std::size_t kk=0;kk<nbOfSubEdges1 && !found;kk++)
 619                     {
 620                       found=(subEdge1PossiblyAlreadyIn1[2*kk]==idBg && subEdge1PossiblyAlreadyIn1[2*kk+1]==idEnd) || (subEdge1PossiblyAlreadyIn1[2*kk]==idEnd && subEdge1PossiblyAlreadyIn1[2*kk+1]==idBg);
 621                       if(!found)
 622                         offset2++;
 623                     }
 624                   if(found)
 625                     {
 626                       ElementaryEdge *e=pol1[offset1+(direct1?offset2:nbOfSubEdges1-offset2-1)];
 627                       e->getPtr()->declareOn();
 628                     }
 629                 }
 630             }
 631           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
 632         }
 633     }
 634 }
 635
 636 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
 637 {
 638   int nbOfNodesInPg=0;
 639   bool presenceOfQuadratic=presenceOfQuadraticEdge();
 640   conn.push_back(presenceOfQuadratic?NORM_QPOLYG:NORM_POLYGON);
 641   for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++)
 642     {
 643       Node *tmp=0;
 644       tmp=(*it)->getStartNode();
 645       std::map<INTERP_KERNEL::Node *,int>::const_iterator it1=mapp.find(tmp);
 646       conn.push_back((*it1).second);
 647       nbOfNodesInPg++;
 648     }
 649   if(presenceOfQuadratic)
 650     {
 651       int j=0;
 652       int off=offset+((int)addCoordsQuadratic.size())/2;
 653       for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++,j++,nbOfNodesInPg++)
 654         {
 655           INTERP_KERNEL::Node *node=(*it)->getPtr()->buildRepresentantOfMySelf();
 656           node->unApplySimilarity(xBary,yBary,fact);
 657           addCoordsQuadratic.push_back((*node)[0]);
 658           addCoordsQuadratic.push_back((*node)[1]);
 659           conn.push_back(off+j);
 660           node->decrRef();
 661         }
 662     }
 663   connI.push_back(connI.back()+nbOfNodesInPg+1);
 664 }
 665
 666 /*!
 667  * 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.
 668  * 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.
 669  * @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
 670  * @param [in,out] edgesBoundaryOther, parameter that stores all edges in result of intersection that are not
 671  */
 672 void QuadraticPolygon::buildPartitionsAbs(QuadraticPolygon& other, std::set<INTERP_KERNEL::Edge *>& edgesThis, std::set<INTERP_KERNEL::Edge *>& edgesBoundaryOther,
 673                                           const std::map<INTERP_KERNEL::Node *,int>& mapp, int idThis, int idOther, int offset,
 674                                           std::vector<double>& addCoordsQuadratic, std::vector<int>& conn, std::vector<int>& connI,
 675                                           std::vector<int>& nbThis, std::vector<int>& nbOther)
 676 {
 677   double xBaryBB, yBaryBB;
 678   double fact=normalizeExt(&other, xBaryBB, yBaryBB);
 679   //Locate \a this relative to \a other (edges of \a this, aka \a pol1 are marked as IN or OUT)
 680   other.performLocatingOperationSlow(*this);  // without any assumption
 681   std::vector<QuadraticPolygon *> res=buildIntersectionPolygons(*this,other);
 682   for(std::vector<QuadraticPolygon *>::iterator it=res.begin();it!=res.end();it++)
 683     {
 684       (*it)->appendCrudeData(mapp,xBaryBB,yBaryBB,fact,offset,addCoordsQuadratic,conn,connI);
 685       INTERP_KERNEL::IteratorOnComposedEdge it1(*it);
 686       for(it1.first();!it1.finished();it1.next())
 687         {
 688           Edge *e=it1.current()->getPtr();
 689           if(edgesThis.find(e)!=edgesThis.end())
 690             edgesThis.erase(e);
 691           else
 692             {
 693               if(edgesBoundaryOther.find(e)!=edgesBoundaryOther.end())
 694                 edgesBoundaryOther.erase(e);
 695               else
 696                 edgesBoundaryOther.insert(e);
 697             }
 698         }
 699       nbThis.push_back(idThis);
 700       nbOther.push_back(idOther);
 701       delete *it;
 702     }
 703   unApplyGlobalSimilarityExt(other,xBaryBB,yBaryBB,fact);
 704 }
 705
 706 /*!
 707  * Remove the two 'identical' edges from the list, and drecRef the objects.
 708  */
 709 void QuadraticPolygon::cleanDegeneratedConsecutiveEdges()
 710 {
 711   IteratorOnComposedEdge it2ii(this);
 712   ElementaryEdge * prevEdge = 0;
 713   int kk = 0;
 714   if  (recursiveSize() > 2)
 715     for(it2ii.first();!it2ii.finished();it2ii.next())
 716       {
 717         ElementaryEdge * cur = it2ii.current();
 718         if (prevEdge && prevEdge->hasSameExtremities(*cur))
 719           {
 720             // Delete the two 'identical' edges:
 721             it2ii.previousLoop(); it2ii.previousLoop();
 722             erase(--kk); erase(kk);
 723             prevEdge = it2ii.current();
 724           }
 725         else
 726           {
 727             kk++;
 728             prevEdge = cur;
 729           }
 730       }
 731 }
 732
 733 /*!
 734  * Warning This method is \b NOT const. 'this' and 'other' are modified after call of this method.
 735  * 'other' is a QuadraticPolygon of \b non closed edges.
 736  */
 737 double QuadraticPolygon::intersectWithAbs1D(QuadraticPolygon& other, bool& isColinear)
 738 {
 739   double ret = 0., xBaryBB, yBaryBB;
 740   double fact = normalize(&other, xBaryBB, yBaryBB);
 741
 742   QuadraticPolygon cpyOfThis(*this);
 743   QuadraticPolygon cpyOfOther(other);
 744   int nbOfSplits = 0;
 745   SplitPolygonsEachOther(cpyOfThis, cpyOfOther, nbOfSplits);
 746   //At this point cpyOfThis and cpyOfOther have been splited at maximum edge so that in/out can been done.
 747   performLocatingOperation(cpyOfOther);
 748   isColinear = false;
 749   for(std::list<ElementaryEdge *>::const_iterator it=cpyOfOther._sub_edges.begin();it!=cpyOfOther._sub_edges.end();it++)
 750     {
 751       switch((*it)->getLoc())
 752       {
 753         case FULL_IN_1:
 754           {
 755             ret += fabs((*it)->getPtr()->getCurveLength());
 756             break;
 757           }
 758         case FULL_ON_1:
 759           {
 760             isColinear=true;
 761             ret += fabs((*it)->getPtr()->getCurveLength());
 762             break;
 763           }
 764         default:
 765           {
 766           }
 767       }
 768     }
 769   return ret * fact;
 770 }
 771
 772 /*!
 773  * Warning contrary to intersectWith method this method is \b NOT const. 'this' and 'other' are modified after call of this method.
 774  */
 775 double QuadraticPolygon::intersectWithAbs(QuadraticPolygon& other, double* barycenter)
 776 {
 777   double ret=0.,bary[2],area,xBaryBB,yBaryBB;
 778   barycenter[0] = barycenter[1] = 0.;
 779   double fact=normalize(&other,xBaryBB,yBaryBB);
 780   std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
 781   for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
 782     {
 783       area=fabs((*iter)->getArea());
 784       (*iter)->getBarycenter(bary);
 785       delete *iter;
 786       ret+=area;
 787       barycenter[0] += bary[0]*area;
 788       barycenter[1] += bary[1]*area;
 789     }
 790   if ( ret > std::numeric_limits<double>::min() )
 791     {
 792       barycenter[0]=barycenter[0]/ret*fact+xBaryBB;
 793       barycenter[1]=barycenter[1]/ret*fact+yBaryBB;
 794
 795     }
 796   return ret*fact*fact;
 797 }
 798
 799 /*!
 800  * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
 801  * This is possible because loc attribute in Edge class is mutable.
 802  * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
 803  */
 804 double QuadraticPolygon::intersectWith(const QuadraticPolygon& other) const
 805 {
 806   double ret=0.;
 807   std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
 808   for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
 809     {
 810       ret+=fabs((*iter)->getArea());
 811       delete *iter;
 812     }
 813   return ret;
 814 }
 815
 816 /*!
 817  * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
 818  * This is possible because loc attribute in Edge class is mutable.
 819  * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
 820  */
 821 double QuadraticPolygon::intersectWith(const QuadraticPolygon& other, double* barycenter) const
 822 {
 823   double ret=0., bary[2];
 824   barycenter[0] = barycenter[1] = 0.;
 825   std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
 826   for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
 827     {
 828       double area = fabs((*iter)->getArea());
 829       (*iter)->getBarycenter(bary);
 830       delete *iter;
 831       ret+=area;
 832       barycenter[0] += bary[0]*area;
 833       barycenter[1] += bary[1]*area;
 834     }
 835   if ( ret > std::numeric_limits<double>::min() )
 836     {
 837       barycenter[0] /= ret;
 838       barycenter[1] /= ret;
 839     }
 840   return ret;
 841 }
 842
 843 /*!
 844  * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
 845  * This is possible because loc attribute in Edge class is mutable.
 846  * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
 847  */
 848 void QuadraticPolygon::intersectForPerimeter(const QuadraticPolygon& other, double& perimeterThisPart, double& perimeterOtherPart, double& perimeterCommonPart) const
 849 {
 850   perimeterThisPart=0.; perimeterOtherPart=0.; perimeterCommonPart=0.;
 851   QuadraticPolygon cpyOfThis(*this);
 852   QuadraticPolygon cpyOfOther(other); int nbOfSplits=0;
 853   SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits);
 854   performLocatingOperation(cpyOfOther);
 855   other.performLocatingOperation(cpyOfThis);
 856   cpyOfThis.dispatchPerimeterExcl(perimeterThisPart,perimeterCommonPart);
 857   cpyOfOther.dispatchPerimeterExcl(perimeterOtherPart,perimeterCommonPart);
 858   perimeterCommonPart/=2.;
 859 }
 860
 861 /*!
 862  * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
 863  * This is possible because loc attribute in Edge class is mutable.
 864  * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
 865  *
 866  * polThis.size()==this->size() and polOther.size()==other.size().
 867  * For each ElementaryEdge of 'this', the corresponding contribution in resulting polygon is in 'polThis'.
 868  * For each ElementaryEdge of 'other', the corresponding contribution in resulting polygon is in 'polOther'.
 869  * As consequence common part are counted twice (in polThis \b and in polOther).
 870  */
 871 void QuadraticPolygon::intersectForPerimeterAdvanced(const QuadraticPolygon& other, std::vector< double >& polThis, std::vector< double >& polOther) const
 872 {
 873   polThis.resize(size());
 874   polOther.resize(other.size());
 875   IteratorOnComposedEdge it1(const_cast<QuadraticPolygon *>(this));
 876   int edgeId=0;
 877   for(it1.first();!it1.finished();it1.next(),edgeId++)
 878     {
 879       ElementaryEdge* curE1=it1.current();
 880       QuadraticPolygon cpyOfOther(other);
 881       QuadraticPolygon tmp;
 882       tmp.pushBack(curE1->clone());
 883       int tmp2;
 884       SplitPolygonsEachOther(tmp,cpyOfOther,tmp2);
 885       other.performLocatingOperation(tmp);
 886       tmp.dispatchPerimeter(polThis[edgeId]);
 887     }
 888   //
 889   IteratorOnComposedEdge it2(const_cast<QuadraticPolygon *>(&other));
 890   edgeId=0;
 891   for(it2.first();!it2.finished();it2.next(),edgeId++)
 892     {
 893       ElementaryEdge* curE2=it2.current();
 894       QuadraticPolygon cpyOfThis(*this);
 895       QuadraticPolygon tmp;
 896       tmp.pushBack(curE2->clone());
 897       int tmp2;
 898       SplitPolygonsEachOther(tmp,cpyOfThis,tmp2);
 899       performLocatingOperation(tmp);
 900       tmp.dispatchPerimeter(polOther[edgeId]);
 901     }
 902 }
 903
 904
 905 /*!
 906  * numberOfCreatedPointsPerEdge is resized to the number of edges of 'this'.
 907  * This method returns in ordered maner the number of newly created points per edge.
 908  * This method performs a split process between 'this' and 'other' that gives the result PThis.
 909  * Then for each edges of 'this' this method counts how many edges in Pthis have the same id.
 910  */
 911 void QuadraticPolygon::intersectForPoint(const QuadraticPolygon& other, std::vector< int >& numberOfCreatedPointsPerEdge) const
 912 {
 913   numberOfCreatedPointsPerEdge.resize(size());
 914   IteratorOnComposedEdge it1(const_cast<QuadraticPolygon *>(this));
 915   int edgeId=0;
 916   for(it1.first();!it1.finished();it1.next(),edgeId++)
 917     {
 918       ElementaryEdge* curE1=it1.current();
 919       QuadraticPolygon cpyOfOther(other);
 920       QuadraticPolygon tmp;
 921       tmp.pushBack(curE1->clone());
 922       int tmp2;
 923       SplitPolygonsEachOther(tmp,cpyOfOther,tmp2);
 924       numberOfCreatedPointsPerEdge[edgeId]=tmp.recursiveSize()-1;
 925     }
 926 }
 927
 928 /*!
 929  * \b WARNING this method is const and other is const too. \b BUT location of Edges in 'this' and 'other' are nevertheless modified.
 930  * This is possible because loc attribute in Edge class is mutable.
 931  * This implies that if 'this' or/and 'other' are reused for intersect* method initLocations has to be called on each of this/them.
 932  * This method is currently not used by any high level functionality.
 933  */
 934 std::vector<QuadraticPolygon *> QuadraticPolygon::intersectMySelfWith(const QuadraticPolygon& other) const
 935 {
 936   QuadraticPolygon cpyOfThis(*this);
 937   QuadraticPolygon cpyOfOther(other); int nbOfSplits=0;
 938   SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits);
 939   //At this point cpyOfThis and cpyOfOther have been splited at maximum edge so that in/out can been done.
 940   performLocatingOperation(cpyOfOther);
 941   return other.buildIntersectionPolygons(cpyOfOther, cpyOfThis);
 942 }
 943
 944 /*!
 945  * This method is typically the first step of boolean operations between pol1 and pol2.
 946  * This method perform the minimal splitting so that at the end each edges constituting pol1 are fully either IN or OUT or ON.
 947  * @param pol1 IN/OUT param that is equal to 'this' when called.
 948  */
 949 void QuadraticPolygon::SplitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticPolygon& pol2, int& nbOfSplits)
 950 {
 951   IteratorOnComposedEdge it1(&pol1),it2(&pol2);
 952   MergePoints merge;
 953   ComposedEdge *c1=new ComposedEdge;
 954   ComposedEdge *c2=new ComposedEdge;
 955   for(it2.first();!it2.finished();it2.next())
 956     {
 957       ElementaryEdge* curE2=it2.current();
 958       if(!curE2->isThereStartPoint())
 959         it1.first();
 960       else
 961         it1=curE2->getIterator();
 962       for(;!it1.finished();)
 963         {
 964
 965           ElementaryEdge* curE1=it1.current();
 966           merge.clear(); nbOfSplits++;
 967           if(curE1->getPtr()->intersectWith(curE2->getPtr(),merge,*c1,*c2))
 968             {
 969               if(!curE1->getDirection()) c1->reverse();
 970               if(!curE2->getDirection()) c2->reverse();
 971               UpdateNeighbours(merge,it1,it2,c1,c2);
 972               //Substitution of simple edge by sub-edges.
 973               delete curE1; // <-- destroying simple edge coming from pol1
 974               delete curE2; // <-- destroying simple edge coming from pol2
 975               it1.insertElemEdges(c1,true);// <-- 2nd param is true to go next.
 976               it2.insertElemEdges(c2,false);// <-- 2nd param is false to avoid to go next.
 977               curE2=it2.current();
 978               //
 979               it1.assignMySelfToAllElems(c2);//To avoid that others
 980               SoftDelete(c1);
 981               SoftDelete(c2);
 982               c1=new ComposedEdge;
 983               c2=new ComposedEdge;
 984             }
 985           else
 986             {
 987               UpdateNeighbours(merge,it1,it2,curE1,curE2);
 988               it1.next();
 989             }
 990         }
 991     }
 992   Delete(c1);
 993   Delete(c2);
 994 }
 995
 996 void QuadraticPolygon::performLocatingOperation(QuadraticPolygon& pol1) const
 997 {
 998   IteratorOnComposedEdge it(&pol1);
 999   TypeOfEdgeLocInPolygon loc=FULL_ON_1;
1000   for(it.first();!it.finished();it.next())
1001     {
1002       ElementaryEdge *cur=it.current();
1003       loc=cur->locateFullyMySelf(*this,loc);//*this=pol2=other
1004     }
1005 }
1006
1007 void QuadraticPolygon::performLocatingOperationSlow(QuadraticPolygon& pol2) const
1008 {
1009   IteratorOnComposedEdge it(&pol2);
1010   for(it.first();!it.finished();it.next())
1011     {
1012       ElementaryEdge *cur=it.current();
1013       cur->locateFullyMySelfAbsolute(*this);
1014     }
1015 }
1016
1017 /*!
1018  * Given 2 polygons \a pol1 and \a pol2 (localized) the resulting polygons are returned.
1019  *
1020  * this : pol1 simplified.
1021  * @param [in] pol1 pol1 split.
1022  * @param [in] pol2 pol2 split.
1023  */
1024 std::vector<QuadraticPolygon *> QuadraticPolygon::buildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2) const
1025 {
1026   std::vector<QuadraticPolygon *> ret;
1027   // Extract from pol1, and pol1 only, all consecutive edges.
1028   // pol1Zip contains concatenated pieces of pol1 which are part of the resulting intersecting cell being built.
1029   std::list<QuadraticPolygon *> pol1Zip=pol1.zipConsecutiveInSegments();
1030   if(!pol1Zip.empty())
1031     ClosePolygons(pol1Zip,*this,pol2,ret);
1032   else
1033     {//borders of pol1 do not cross pol2,and pol1 borders are outside of pol2. That is to say, either pol1 and pol2
1034       //do not overlap or  pol2 is fully inside pol1. So in the first case no intersection, in the other case
1035       //the intersection is pol2.
1036       ElementaryEdge *e1FromPol2=pol2[0];
1037       TypeOfEdgeLocInPolygon loc=FULL_ON_1;
1038       loc=e1FromPol2->locateFullyMySelf(*this,loc);
1039       if(loc==FULL_IN_1)
1040         ret.push_back(new QuadraticPolygon(pol2));
1041     }
1042   return ret;
1043 }
1044
1045 /*!
1046  * Returns parts of potentially non closed-polygons. Each returned polygons are not mergeable.
1047  * this : pol1 split and localized.
1048  */
1049 std::list<QuadraticPolygon *> QuadraticPolygon::zipConsecutiveInSegments() const
1050 {
1051   std::list<QuadraticPolygon *> ret;
1052   IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(this));
1053   int nbOfTurns=recursiveSize();
1054   int i=0;
1055   if(!it.goToNextInOn(false,i,nbOfTurns))
1056     return ret;
1057   i=0;
1058   //
1059   while(i<nbOfTurns)
1060     {
1061       QuadraticPolygon *tmp1=new QuadraticPolygon;
1062       TypeOfEdgeLocInPolygon loc=it.current()->getLoc();
1063       while(loc!=FULL_OUT_1 && i<nbOfTurns)
1064         {
1065           ElementaryEdge *tmp3=it.current()->clone();
1066           tmp1->pushBack(tmp3);
1067           it.nextLoop(); i++;
1068           loc=it.current()->getLoc();
1069         }
1070       if(tmp1->empty())
1071         {
1072           delete tmp1;
1073           continue;
1074         }
1075       ret.push_back(tmp1);
1076       it.goToNextInOn(true,i,nbOfTurns);
1077     }
1078   return ret;
1079 }
1080
1081 /*!
1082  * @param [in] pol1zip is a list of set of edges (=an opened polygon) coming from split polygon 1.
1083  * @param [in] pol1 should be considered as pol1Simplified.
1084  * @param [in] pol2 is split pol2.
1085  * @param [out] results the resulting \b CLOSED polygons.
1086  */
1087 void QuadraticPolygon::ClosePolygons(std::list<QuadraticPolygon *>& pol1Zip, const QuadraticPolygon& pol1, const QuadraticPolygon& pol2,
1088                                      std::vector<QuadraticPolygon *>& results)
1089 {
1090   bool directionKnownInPol2=false;
1091   bool directionInPol2;
1092   for(std::list<QuadraticPolygon *>::iterator iter=pol1Zip.begin();iter!=pol1Zip.end();)
1093     {
1094       // Build incrementally the full closed cells from the consecutive line parts already built in pol1Zip.
1095       // At the end of the process the item currently iterated has been totally completed (start_node=end_node)
1096       // This process can produce several cells.
1097       if((*iter)->completed())
1098         {
1099           results.push_back(*iter);
1100           directionKnownInPol2=false;
1101           iter=pol1Zip.erase(iter);
1102           continue;
1103         }
1104       if(!directionKnownInPol2)
1105         {
1106           if(!(*iter)->haveIAChanceToBeCompletedBy(pol1,pol2,directionInPol2))
1107             { delete *iter; iter=pol1Zip.erase(iter); continue; }
1108           else
1109             directionKnownInPol2=true;
1110         }
1111       std::list<QuadraticPolygon *>::iterator iter2=iter; iter2++;
1112       // Fill as much as possible the current iterate (=a part of pol1) with consecutive pieces from pol2:
1113       std::list<QuadraticPolygon *>::iterator iter3=(*iter)->fillAsMuchAsPossibleWith(pol2,iter2,pol1Zip.end(),directionInPol2);
1114       // and now add a full connected piece from pol1Zip:
1115       if(iter3!=pol1Zip.end())
1116         {
1117           (*iter)->pushBack(*iter3);
1118           SoftDelete(*iter3);
1119           pol1Zip.erase(iter3);
1120         }
1121     }
1122 }
1123
1124 /*!
1125  * 'this' is expected to be set of edges (not closed) of pol1 split.
1126  */
1127 bool QuadraticPolygon::haveIAChanceToBeCompletedBy(const QuadraticPolygon& pol1NotSplitted, const QuadraticPolygon& pol2Splitted, bool& direction) const
1128 {
1129   IteratorOnComposedEdge it2(const_cast<QuadraticPolygon *>(&pol2Splitted));
1130   bool found=false;
1131   Node *n=getEndNode();
1132   ElementaryEdge *cur=it2.current();
1133   for(it2.first();!it2.finished() && !found;)
1134     {
1135       cur=it2.current();
1136       found=(cur->getStartNode()==n);
1137       if(!found)
1138         it2.next();
1139     }
1140   if(!found)
1141     throw Exception("Internal error: polygons incompatible with each others. Should never happen!");
1142   //Ok we found correspondence between this and pol2. Searching for right direction to close polygon.
1143   ElementaryEdge *e=_sub_edges.back();
1144   if(e->getLoc()==FULL_ON_1)
1145     {
1146       if(e->getPtr()==cur->getPtr())
1147         {
1148           direction=false;
1149           it2.previousLoop();
1150           cur=it2.current();
1151           Node *repr=cur->getPtr()->buildRepresentantOfMySelf();
1152           bool ret=pol1NotSplitted.isInOrOut(repr);
1153           repr->decrRef();
1154           return ret;
1155         }
1156       else
1157         {
1158           direction=true;
1159           Node *repr=cur->getPtr()->buildRepresentantOfMySelf();
1160           bool ret=pol1NotSplitted.isInOrOut(repr);
1161           repr->decrRef();
1162           return ret;
1163         }
1164     }
1165   else
1166     direction=cur->locateFullyMySelfAbsolute(pol1NotSplitted)==FULL_IN_1;
1167   return true;
1168 }
1169
1170 /*!
1171  * This method fills as much as possible \a this (a sub-part of pol1 split) with edges of \a pol2Splitted.
1172  */
1173 std::list<QuadraticPolygon *>::iterator QuadraticPolygon::fillAsMuchAsPossibleWith(const QuadraticPolygon& pol2Splitted,
1174                                                                                    std::list<QuadraticPolygon *>::iterator iStart,
1175                                                                                    std::list<QuadraticPolygon *>::iterator iEnd,
1176                                                                                    bool direction)
1177 {
1178   IteratorOnComposedEdge it1(const_cast<QuadraticPolygon *>(&pol2Splitted));
1179   bool found=false;
1180   Node *n=getEndNode();
1181   ElementaryEdge *cur1;
1182   for(it1.first();!it1.finished() && !found;)
1183     {
1184       cur1=it1.current();
1185       found=(cur1->getStartNode()==n);
1186       if(!found)
1187         it1.next();
1188     }
1189   if(!direction)
1190     it1.previousLoop();
1191   Node *nodeToTest;
1192   int szMax(pol2Splitted.size()+1),ii(0);   // protection against aggressive users of IntersectMeshes using invalid input meshes ...
1193   std::list<QuadraticPolygon *>::iterator ret;
1194   do
1195     { // Stack (consecutive) edges of pol1 into the result (no need to care about ordering, edges from pol1 are already consecutive)
1196       cur1=it1.current();
1197       ElementaryEdge *tmp=cur1->clone();
1198       if(!direction)
1199         tmp->reverse();
1200       pushBack(tmp);
1201       nodeToTest=tmp->getEndNode();
1202       direction?it1.nextLoop():it1.previousLoop();
1203       ret=CheckInList(nodeToTest,iStart,iEnd);
1204       if(completed())
1205         return iEnd;
1206       ii++;
1207     }
1208   while(ret==iEnd && ii<szMax);
1209   if(ii==szMax)// here a protection against aggressive users of IntersectMeshes of invalid input meshes
1210     throw INTERP_KERNEL::Exception("QuadraticPolygon::fillAsMuchAsPossibleWith : Something is invalid with input polygons !");
1211   return ret;
1212 }
1213
1214 std::list<QuadraticPolygon *>::iterator QuadraticPolygon::CheckInList(Node *n, std::list<QuadraticPolygon *>::iterator iStart,
1215                                                                       std::list<QuadraticPolygon *>::iterator iEnd)
1216 {
1217   for(std::list<QuadraticPolygon *>::iterator iter=iStart;iter!=iEnd;iter++)
1218     if((*iter)->isNodeIn(n))
1219       return iter;
1220   return iEnd;
1221 }
1222
1223 /*!
1224 * Compute the remaining parts of the intersection of mesh1 by mesh2.
1225 * The general algorithm is to :
1226 *   - either return full cells from pol1 that were simply not touched by mesh2
1227 *   - or to:
1228 *      - concatenate pieces from pol1 into consecutive pieces (a bit like zipConsecutiveSegments())
1229 *      - complete those pieces by edges found in edgesInPol2OnBoundary, which are edges from pol2 located on the boundary of the previously built
1230 *      intersecting cells
1231 */
1232 void QuadraticPolygon::ComputeResidual(const QuadraticPolygon& pol1, const std::set<Edge *>& notUsedInPol1, const std::set<Edge *>& edgesInPol2OnBoundary,
1233                                        const std::map<INTERP_KERNEL::Node *,int>& mapp, int offset, int idThis,
1234                                        std::vector<double>& addCoordsQuadratic, std::vector<int>& conn,
1235                                        std::vector<int>& connI, std::vector<int>& nb1, std::vector<int>& nb2)
1236 {
1237   // Initialise locations on pol1. Remember that edges found in 'notUsedInPol1' are also part of the edges forming pol1.
1238   pol1.initLocations();
1239   for(std::set<Edge *>::const_iterator it1=notUsedInPol1.begin();it1!=notUsedInPol1.end();it1++)
1240     { (*it1)->initLocs(); (*it1)->declareOn(); }
1241   for(std::set<Edge *>::const_iterator it2=edgesInPol2OnBoundary.begin();it2!=edgesInPol2OnBoundary.end();it2++)
1242     { (*it2)->initLocs(); (*it2)->declareIn(); }
1243   ////
1244   std::set<Edge *> notUsedInPol1L(notUsedInPol1);
1245   IteratorOnComposedEdge itPol1(const_cast<QuadraticPolygon *>(&pol1));
1246   int sz=pol1.size();
1247   std::list<QuadraticPolygon *> pol1Zip;
1248   // If none of the edges of pol1 was consumed by the rebuilding process, we can directly take pol1 as it is to form a cell:
1249   if(pol1.size()==(int)notUsedInPol1.size() && edgesInPol2OnBoundary.empty())
1250     {
1251       pol1.appendCrudeData(mapp,0.,0.,1.,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1);
1252       return ;
1253     }
1254   // Zip consecutive edges found in notUsedInPol1L and which are not overlapping with boundary edge from edgesInPol2OnBoundary:
1255   while(!notUsedInPol1L.empty())
1256     {
1257       // If all nodes are IN or edge is ON (i.e. as initialised at the begining of the method) then error
1258       for(int i=0;i<sz && (itPol1.current()->getStartNode()->getLoc()!=IN_1 || itPol1.current()->getLoc()!=FULL_ON_1);i++)
1259         itPol1.nextLoop();
1260       if(itPol1.current()->getStartNode()->getLoc()!=IN_1 || itPol1.current()->getLoc()!=FULL_ON_1)
1261         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 !");
1262       QuadraticPolygon *tmp1=new QuadraticPolygon;
1263       do
1264         {
1265           Edge *ee=itPol1.current()->getPtr();
1266           if(ee->getLoc()==FULL_ON_1)
1267             {
1268               ee->incrRef(); notUsedInPol1L.erase(ee);
1269               tmp1->pushBack(new ElementaryEdge(ee,itPol1.current()->getDirection()));
1270             }
1271           itPol1.nextLoop();
1272         }
1273       while(itPol1.current()->getStartNode()->getLoc()!=IN_1 && !notUsedInPol1L.empty());
1274       pol1Zip.push_back(tmp1);
1275     }
1276
1277   ////
1278   std::list<QuadraticPolygon *> retPolsUnderContruction;
1279   std::list<Edge *> edgesInPol2OnBoundaryL(edgesInPol2OnBoundary.begin(),edgesInPol2OnBoundary.end());
1280   std::map<QuadraticPolygon *, std::list<QuadraticPolygon *> > pol1ZipConsumed;  // for memory management only.
1281   std::size_t maxNbOfTurn=edgesInPol2OnBoundaryL.size(),nbOfTurn=0,iiMNT=0;
1282   for(std::list<QuadraticPolygon *>::const_iterator itMNT=pol1Zip.begin();itMNT!=pol1Zip.end();itMNT++,iiMNT++)
1283     nbOfTurn+=(*itMNT)->size();
1284   maxNbOfTurn=maxNbOfTurn*nbOfTurn; maxNbOfTurn*=maxNbOfTurn;
1285   // [ABN] at least 3 turns for very small cases (typically one (quad) edge against one (quad or lin) edge forming a new cell)!
1286   maxNbOfTurn = maxNbOfTurn<3 ? 3 : maxNbOfTurn;
1287   nbOfTurn=0;
1288   while(nbOfTurn<maxNbOfTurn)  // the 'normal' way out of this loop is the break towards the end when pol1Zip is empty.
1289     {
1290       // retPolsUnderConstruction initially empty -> see if(!pol1Zip.empty()) below ...
1291       for(std::list<QuadraticPolygon *>::iterator itConstr=retPolsUnderContruction.begin();itConstr!=retPolsUnderContruction.end();)
1292         {
1293           if((*itConstr)->getStartNode()==(*itConstr)->getEndNode())  // reconstruction of a cell is finished
1294             {
1295               itConstr++;
1296               continue;
1297             }
1298           Node *curN=(*itConstr)->getEndNode();
1299           bool smthHappened=false;
1300           // Complete a partially reconstructed polygon with boundary edges by matching nodes:
1301           for(std::list<Edge *>::iterator it2=edgesInPol2OnBoundaryL.begin();it2!=edgesInPol2OnBoundaryL.end();)
1302             {
1303               if(curN==(*it2)->getStartNode())
1304                 { (*it2)->incrRef(); (*itConstr)->pushBack(new ElementaryEdge(*it2,true)); curN=(*it2)->getEndNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); }
1305               else if(curN==(*it2)->getEndNode())
1306                 { (*it2)->incrRef(); (*itConstr)->pushBack(new ElementaryEdge(*it2,false)); curN=(*it2)->getStartNode(); smthHappened=true; it2=edgesInPol2OnBoundaryL.erase(it2); }
1307               else
1308                 it2++;
1309             }
1310           if(smthHappened)
1311             {
1312               for(std::list<QuadraticPolygon *>::iterator itZip=pol1Zip.begin();itZip!=pol1Zip.end();)
1313                 {
1314                   if(curN==(*itZip)->getStartNode()) // we found a matching piece to append in pol1Zip. Append all of it to the current polygon being reconstr
1315                     {
1316                       for(std::list<ElementaryEdge *>::const_iterator it4=(*itZip)->_sub_edges.begin();it4!=(*itZip)->_sub_edges.end();it4++)
1317                         { (*it4)->getPtr()->incrRef(); bool dir=(*it4)->getDirection(); (*itConstr)->pushBack(new ElementaryEdge((*it4)->getPtr(),dir)); }
1318                       pol1ZipConsumed[*itConstr].push_back(*itZip);
1319                       curN=(*itZip)->getEndNode();
1320                       itZip=pol1Zip.erase(itZip);  // one zipped piece has been consumed
1321                       break;                       // we can stop here, pieces in pol1Zip are not connected, by definition.
1322                     }
1323                   else
1324                     itZip++;
1325                 }
1326             }
1327           else
1328             {
1329               for(std::list<ElementaryEdge *>::const_iterator it5=(*itConstr)->_sub_edges.begin();it5!=(*itConstr)->_sub_edges.end();it5++)
1330                 {
1331                   Edge *ee=(*it5)->getPtr();
1332                   if(edgesInPol2OnBoundary.find(ee)!=edgesInPol2OnBoundary.end())
1333                     edgesInPol2OnBoundaryL.push_back(ee);
1334                 }
1335               for(std::list<QuadraticPolygon *>::iterator it6=pol1ZipConsumed[*itConstr].begin();it6!=pol1ZipConsumed[*itConstr].end();it6++)
1336                 pol1Zip.push_front(*it6);
1337               pol1ZipConsumed.erase(*itConstr);
1338               delete *itConstr;
1339               itConstr=retPolsUnderContruction.erase(itConstr);
1340             }
1341         }
1342       if(!pol1Zip.empty())  // the filling process of retPolsUnderConstruction starts here
1343         {
1344           QuadraticPolygon *tmp=new QuadraticPolygon;
1345           QuadraticPolygon *first=*(pol1Zip.begin());
1346           for(std::list<ElementaryEdge *>::const_iterator it4=first->_sub_edges.begin();it4!=first->_sub_edges.end();it4++)
1347             { (*it4)->getPtr()->incrRef(); bool dir=(*it4)->getDirection(); tmp->pushBack(new ElementaryEdge((*it4)->getPtr(),dir)); }
1348           pol1ZipConsumed[tmp].push_back(first);
1349           retPolsUnderContruction.push_back(tmp);
1350           pol1Zip.erase(pol1Zip.begin());
1351         }
1352       else
1353         break;
1354       nbOfTurn++;
1355     }
1356   if(nbOfTurn==maxNbOfTurn)
1357     {
1358       std::ostringstream oss; oss << "Error during reconstruction of residual of cell ! It appears that either source or/and target mesh is/are not conform !";
1359       oss << " Number of turns is = " << nbOfTurn << " !";
1360       throw INTERP_KERNEL::Exception(oss.str().c_str());
1361     }
1362   // Convert to integer connectivity:
1363   for(std::list<QuadraticPolygon *>::iterator itConstr=retPolsUnderContruction.begin();itConstr!=retPolsUnderContruction.end();itConstr++)
1364     {
1365       if((*itConstr)->getStartNode()==(*itConstr)->getEndNode())  // take only fully closed reconstructed polygon (?? might there be others??)
1366         {
1367           (*itConstr)->appendCrudeData(mapp,0.,0.,1.,offset,addCoordsQuadratic,conn,connI); nb1.push_back(idThis); nb2.push_back(-1);
1368           for(std::list<QuadraticPolygon *>::iterator it6=pol1ZipConsumed[*itConstr].begin();it6!=pol1ZipConsumed[*itConstr].end();it6++)
1369             delete *it6;
1370           delete *itConstr;
1371         }
1372       else
1373         {
1374           std::ostringstream oss; oss << "Internal error during reconstruction of residual of cell! Non fully closed polygon built!";
1375           throw INTERP_KERNEL::Exception(oss.str().c_str());
1376         }
1377     }
1378 }