src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx

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