src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DEdge.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 "InterpKernelGeo2DEdge.hxx"
  22 #include "InterpKernelGeo2DEdgeLin.hxx"
  23 #include "InterpKernelGeo2DEdgeInfLin.hxx"
  24 //#include "EdgeParabol.hxx"
  25 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
  26 #include "InterpKernelGeo2DComposedEdge.hxx"
  27 #include "InterpKernelException.hxx"
  28
  29 #include <algorithm>
  30
  31 #include <functional>
  32
  33
  34 using namespace INTERP_KERNEL;
  35
  36 MergePoints::MergePoints():_ass1Start1(0),_ass1End1(0),_ass1Start2(0),_ass1End2(0),
  37     _ass2Start1(0),_ass2End1(0),_ass2Start2(0),_ass2End2(0)
  38 {
  39 }
  40
  41 void MergePoints::start1Replaced()
  42 {
  43   unsigned nbOfAsso=getNumberOfAssociations();
  44   if(nbOfAsso==0)
  45     _ass1Start1=1;
  46   else
  47     _ass2Start1=1;
  48 }
  49
  50 void MergePoints::end1Replaced()
  51 {
  52   unsigned nbOfAsso=getNumberOfAssociations();
  53   if(nbOfAsso==0)
  54     _ass1End1=1;
  55   else
  56     _ass2End1=1;
  57 }
  58
  59 void MergePoints::start1OnStart2()
  60 {
  61   unsigned nbOfAsso=getNumberOfAssociations();
  62   if(nbOfAsso==0)
  63     {
  64       _ass1Start1=1;
  65       _ass1Start2=1;
  66     }
  67   else
  68     {
  69       _ass2Start1=1;
  70       _ass2Start2=1;
  71     }
  72 }
  73
  74 void MergePoints::start1OnEnd2()
  75 {
  76   unsigned nbOfAsso=getNumberOfAssociations();
  77   if(nbOfAsso==0)
  78     {
  79       _ass1Start1=1;
  80       _ass1End2=1;
  81     }
  82   else
  83     {
  84       _ass2Start1=1;
  85       _ass2End2=1;
  86     }
  87 }
  88
  89 void MergePoints::end1OnStart2()
  90 {
  91   unsigned nbOfAsso=getNumberOfAssociations();
  92   if(nbOfAsso==0)
  93     {
  94       _ass1End1=1;
  95       _ass1Start2=1;
  96     }
  97   else
  98     {
  99       _ass2End1=1;
 100       _ass2Start2=1;
 101     }
 102 }
 103
 104 void MergePoints::end1OnEnd2()
 105 {
 106   unsigned nbOfAsso=getNumberOfAssociations();
 107   if(nbOfAsso==0)
 108     {
 109       _ass1End1=1;
 110       _ass1End2=1;
 111     }
 112   else
 113     {
 114       _ass2End1=1;
 115       _ass2End2=1;
 116     }
 117 }
 118
 119 bool MergePoints::isStart1(unsigned rk) const
 120 {
 121   if(rk==0)
 122     return _ass1Start1;
 123   else
 124     return _ass2Start1;
 125 }
 126
 127 bool MergePoints::isEnd1(unsigned rk) const
 128 {
 129   if(rk==0)
 130     return _ass1End1;
 131   else
 132     return _ass2End1;
 133 }
 134
 135 bool MergePoints::isStart2(unsigned rk) const
 136 {
 137   if(rk==0)
 138     return _ass1Start2;
 139   else
 140     return _ass2Start2;
 141 }
 142
 143 bool MergePoints::isEnd2(unsigned rk) const
 144 {
 145   if(rk==0)
 146     return _ass1End2;
 147   else
 148     return _ass2End2;
 149 }
 150
 151 void MergePoints::clear()
 152 {
 153   _ass1Start1=0;_ass1End1=0;_ass1Start2=0;_ass1End2=0;
 154   _ass2Start1=0;_ass2End1=0;_ass2Start2=0;_ass2End2=0;
 155 }
 156
 157 unsigned MergePoints::getNumberOfAssociations() const
 158 {
 159   unsigned ret=0;
 160   unsigned subTot=_ass1Start1+_ass1End1+_ass1Start2+_ass1End2;
 161   if(subTot!=0)
 162     ret++;
 163   subTot=_ass2Start1+_ass2End1+_ass2Start2+_ass2End2;
 164   if(subTot!=0)
 165     ret++;
 166   return ret;
 167 }
 168
 169 void MergePoints::PushInMap(int key, int value, std::map<int,int>& mergedNodes)
 170 {
 171   if(key!=-1 && value!=-1)
 172     mergedNodes[key]=value;
 173 }
 174
 175 void MergePoints::updateMergedNodes(int e1Start, int e1End, int e2Start, int e2End, std::map<int,int>& mergedNodes)
 176 {
 177   unsigned subTot(_ass1Start1+_ass1End1+_ass1Start2+_ass1End2);
 178   if(subTot!=0)
 179     {
 180       if(_ass1Start1 && _ass1Start2)
 181         PushInMap(e2Start,e1Start,mergedNodes);
 182       if(_ass1Start1 && _ass1End2)
 183         PushInMap(e2End,e1Start,mergedNodes);
 184       if(_ass1End1 && _ass1Start2)
 185         PushInMap(e2Start,e1End,mergedNodes);
 186       if(_ass1End1 && _ass1End2)
 187         PushInMap(e2End,e1End,mergedNodes);
 188     }
 189   subTot=_ass2Start1+_ass2End1+_ass2Start2+_ass2End2;
 190   if(subTot!=0)
 191     {
 192       if(_ass2Start1 && _ass2Start2)
 193         PushInMap(e2Start,e1Start,mergedNodes);
 194       if(_ass2Start1 && _ass2End2)
 195         PushInMap(e2End,e1Start,mergedNodes);
 196       if(_ass2End1 && _ass2Start2)
 197         PushInMap(e2Start,e1End,mergedNodes);
 198       if(_ass2End1 && _ass2End2)
 199         PushInMap(e2End,e1End,mergedNodes);
 200     }
 201 }
 202
 203 IntersectElement::IntersectElement(double val1, double val2, bool start1, bool end1, bool start2, bool end2, Node *node
 204                                    , const Edge& e1, const Edge& e2, bool keepOrder):_1S(keepOrder?start1:start2),
 205                                        _1E(keepOrder?end1:end2),
 206                                        _2S(keepOrder?start2:start1),
 207                                        _2E(keepOrder?end2:end1),
 208                                        _chararct_val_for_e1(keepOrder?val1:val2),
 209                                        _chararct_val_for_e2(keepOrder?val2:val1),
 210                                        _node(node),_loc_of_node(node->getLoc()),_e1(keepOrder?e1:e2),
 211                                        _e2(keepOrder?e2:e1)
 212 {
 213 }
 214
 215 IntersectElement::IntersectElement(const IntersectElement& other):_1S(other._1S),_1E(other._1E),_2S(other._2S),_2E(other._2E),
 216     _chararct_val_for_e1(other._chararct_val_for_e1),
 217     _chararct_val_for_e2(other._chararct_val_for_e2),_node(other._node),
 218     _loc_of_node(other._loc_of_node),_e1(other._e1), _e2(other._e2)
 219 {
 220   if(_node)
 221     _node->incrRef();
 222 }
 223
 224 IntersectElement& IntersectElement::operator=(const IntersectElement& other)
 225 {
 226   _1S=other._1S;_1E=other._1E; _2S=other._2S; _2E=other._2E;
 227   _chararct_val_for_e1=other._chararct_val_for_e1;
 228   _chararct_val_for_e2=other._chararct_val_for_e2;
 229   setNode(other._node);
 230   return *this;
 231 }
 232
 233 bool IntersectElement::operator<(const IntersectElement& other) const
 234 {
 235   return _e1.isLower(_chararct_val_for_e1,other._chararct_val_for_e1);
 236 }
 237
 238 IntersectElement::~IntersectElement()
 239 {
 240   if(_node)
 241     _node->decrRef();
 242 }
 243
 244 /*!
 245  * Returns 0 or 1.
 246  */
 247 bool IntersectElement::isOnMergedExtremity() const
 248 {
 249   if( (_1S && _2S) || (_1S && _2E) || (_1E && _2S) || (_1E && _2E) )
 250     return true;
 251   return false;
 252 }
 253
 254 /*!
 255  * To call if isOnMergedExtremity returned true.
 256  */
 257 void IntersectElement::performMerging(MergePoints& commonNode) const
 258 {
 259   if(_1S && _2S)
 260     {
 261       if(_e1.changeStartNodeWith(_e2.getStartNode()))
 262         {
 263           _e2.getStartNode()->declareOnLim();
 264           commonNode.start1OnStart2();
 265         }
 266     }
 267   else if(_1S && _2E)
 268     {
 269       if(_e1.changeStartNodeWith(_e2.getEndNode()))
 270         {
 271           _e2.getEndNode()->declareOnLim();
 272           commonNode.start1OnEnd2();
 273         }
 274     }
 275   else if(_1E && _2S)
 276     {
 277       if(_e1.changeEndNodeWith(_e2.getStartNode()))
 278         {
 279           _e2.getStartNode()->declareOnLim();
 280           commonNode.end1OnStart2();
 281         }
 282     }
 283   else if(_1E && _2E)
 284     {
 285       if(_e1.changeEndNodeWith(_e2.getEndNode()))
 286         {
 287           _e2.getEndNode()->declareOnLim();
 288           commonNode.end1OnEnd2();
 289         }
 290     }
 291 }
 292
 293 /*!
 294  * This method is const because 'node' is supposed to be equal geometrically to _node.
 295  */
 296 void IntersectElement::setNode(Node *node) const
 297 {
 298   if(node!=_node)
 299     {
 300       if(_node)
 301         ((Node *)_node)->decrRef();
 302       (const_cast<IntersectElement *>(this))->_node=node;
 303       if(_node)
 304         _node->incrRef();
 305     }
 306 }
 307
 308 bool IntersectElement::isLowerOnOther(const IntersectElement& other) const
 309 {
 310   return _e2.isLower(_chararct_val_for_e2,other._chararct_val_for_e2);
 311 }
 312
 313 unsigned IntersectElement::isOnExtrForAnEdgeAndInForOtherEdge() const
 314 {
 315   if(( _1S && !(_2S || _2E) ) || ( _1E && !(_2S || _2E) ))
 316     {
 317       if(_1S && !(_2S || _2E))
 318         setNode(_e1.getStartNode());
 319       else
 320         setNode(_e1.getEndNode());
 321       if(_e2.isIn(_chararct_val_for_e2))
 322         return LIMIT_ON;
 323       return LIMIT_ALONE;
 324     }
 325   if(( _2S && !(_1S || _1E) ) || ( _2E && !(_1S || _1E)))
 326     {
 327       if(_2S && !(_1S || _1E))
 328         setNode(_e2.getStartNode());
 329       else
 330         setNode(_e2.getEndNode());
 331       if(_e1.isIn(_chararct_val_for_e1))
 332         return LIMIT_ON;
 333       return LIMIT_ALONE;
 334     }
 335   return NO_LIMIT;
 336 }
 337
 338 bool IntersectElement::isIncludedByBoth() const
 339 {
 340   return _e1.isIn(_chararct_val_for_e1) && _e2.isIn(_chararct_val_for_e2);
 341 }
 342
 343 bool EdgeIntersector::intersect(const Bounds *whereToFind, std::vector<Node *>& newNodes, bool& order, MergePoints& commonNode)
 344 {
 345   std::list< IntersectElement > listOfIntesc=getIntersectionsCharacteristicVal();
 346   std::list< IntersectElement >::iterator iter;
 347   for(iter=listOfIntesc.begin();iter!=listOfIntesc.end();)
 348     {
 349       if((*iter).isOnMergedExtremity())
 350         {
 351           (*iter).performMerging(commonNode);
 352           iter=listOfIntesc.erase(iter);
 353           continue;
 354         }
 355       unsigned tmp=(*iter).isOnExtrForAnEdgeAndInForOtherEdge();
 356       if(tmp==IntersectElement::LIMIT_ALONE)
 357         {
 358           iter=listOfIntesc.erase(iter);
 359           continue;
 360         }
 361       else if(tmp==IntersectElement::LIMIT_ON)
 362         {
 363           (*iter).attachLoc();
 364           iter++;
 365           continue;
 366         }
 367       if(!(*iter).isIncludedByBoth())
 368         {
 369           iter=listOfIntesc.erase(iter);
 370           continue;
 371         }
 372       (*iter).attachLoc();
 373       iter++;
 374     }
 375   if(listOfIntesc.size()==0)
 376     return false;
 377   if(listOfIntesc.size()==1)
 378     {
 379       order=true;//useless
 380       newNodes.push_back(listOfIntesc.front().getNodeAndReleaseIt());
 381     }
 382   else
 383     {
 384       std::vector<IntersectElement> vecOfIntesc(listOfIntesc.begin(),listOfIntesc.end());
 385       listOfIntesc.clear();
 386       sort(vecOfIntesc.begin(),vecOfIntesc.end());
 387       for(std::vector<IntersectElement>::iterator iterV=vecOfIntesc.begin();iterV!=vecOfIntesc.end();iterV++)
 388         newNodes.push_back((*iterV).getNodeAndReleaseIt());
 389       order=vecOfIntesc.front().isLowerOnOther(vecOfIntesc.back());
 390     }
 391   return true;
 392 }
 393
 394 /*!
 395  * Locates 'node' regarding edge this->_e1. If node is located close to (with distant lt epsilon) start or end point of _e1,
 396  * 'node' takes its place. In this case 'obvious' is set to true and 'commonNode' stores information of merge point and finally 'where' is set.
 397  * Furthermore 'node' is declared as ON LIMIT to indicate in locating process that an absolute location computation will have to be done.
 398  * If 'node' is not close to start or end point of _e1, 'obvious' is set to false and 'commonNode' and 'where' are let unchanged.
 399  */
 400 void EdgeIntersector::obviousCaseForCurvAbscisse(Node *node, TypeOfLocInEdge& where, MergePoints& commonNode, bool& obvious) const
 401 {
 402   obvious=true;
 403   if(node->isEqual(*_e1.getStartNode()))
 404     {
 405       where=START;
 406       if(_e1.changeStartNodeWith(node))
 407         {
 408           commonNode.start1Replaced();
 409           node->declareOnLim();
 410         }
 411       return ;
 412     }
 413   if(node->isEqual(*_e1.getEndNode()))
 414     {
 415       where=END;
 416       if(_e1.changeEndNodeWith(node))
 417         {
 418           commonNode.end1Replaced();
 419           node->declareOnLim();
 420         }
 421       return ;
 422     }
 423   obvious=false;
 424 }
 425
 426 Edge::Edge(double sX, double sY, double eX, double eY):_cnt(1),_loc(FULL_UNKNOWN),_start(new Node(sX,sY)),_end(new Node(eX,eY))
 427 {
 428 }
 429
 430 Edge::~Edge()
 431 {
 432   _start->decrRef();
 433   if(_end)
 434     _end->decrRef();
 435 }
 436
 437 bool Edge::decrRef()
 438 {
 439   bool ret=(--_cnt==0);
 440   if(ret)
 441     delete this;
 442   return ret;
 443 }
 444
 445 void Edge::declareOn() const
 446 {
 447   if(_loc==FULL_UNKNOWN)
 448     {
 449       _loc=FULL_ON_1;
 450       _start->declareOn();
 451       _end->declareOn();
 452     }
 453 }
 454
 455 void Edge::declareIn() const
 456 {
 457   if(_loc==FULL_UNKNOWN)
 458     {
 459       _loc=FULL_IN_1;
 460       _start->declareIn();
 461       _end->declareIn();
 462     }
 463 }
 464
 465 void Edge::declareOut() const
 466 {
 467   if(_loc==FULL_UNKNOWN)
 468     {
 469       _loc=FULL_OUT_1;
 470       _start->declareOut();
 471       _end->declareOut();
 472     }
 473 }
 474
 475 void Edge::fillXfigStreamForLoc(std::ostream& stream) const
 476 {
 477   switch(_loc)
 478   {
 479     case FULL_IN_1:
 480       stream << '2';//Green
 481       break;
 482     case FULL_OUT_1:
 483       stream << '1';//Bleue
 484       break;
 485     case FULL_ON_1:
 486       stream << '4';//Red
 487       break;
 488     default:
 489       stream << '0';
 490   }
 491 }
 492
 493 bool Edge::changeStartNodeWith(Node *otherStartNode) const
 494 {
 495   if(_start==otherStartNode)
 496     return true;
 497   if(_start->isEqual(*otherStartNode))
 498     {
 499       ((const_cast<Edge *>(this))->_start)->decrRef();//un-const cast Ok thanks to 2 lines above.
 500       ((const_cast<Edge *>(this))->_start)=otherStartNode;
 501       _start->incrRef();
 502       return true;
 503     }
 504   return false;
 505 }
 506
 507 bool Edge::changeStartNodeWithAndKeepTrack(Node *otherStartNode, std::vector<Node *>& track) const
 508 {
 509   if(_start==otherStartNode)
 510     return true;
 511   if(_start->isEqualAndKeepTrack(*otherStartNode,track))
 512     {
 513       ((const_cast<Edge *>(this))->_start)->decrRef();//un-const cast Ok thanks to 2 lines above.
 514       ((const_cast<Edge *>(this))->_start)=otherStartNode;
 515       otherStartNode->incrRef();
 516       return true;
 517     }
 518   return false;
 519 }
 520
 521 bool Edge::changeEndNodeWith(Node *otherEndNode) const
 522 {
 523   if(_end==otherEndNode)
 524     return true;
 525   if(_end->isEqual(*otherEndNode))
 526     {
 527       ((const_cast<Edge *>(this))->_end)->decrRef();
 528       ((const_cast<Edge *>(this))->_end)=otherEndNode;
 529       _end->incrRef();
 530       return true;
 531     }
 532   return false;
 533 }
 534
 535 bool Edge::changeEndNodeWithAndKeepTrack(Node *otherEndNode, std::vector<Node *>& track) const
 536 {
 537   if(_end==otherEndNode)
 538     return true;
 539   if(_end->isEqualAndKeepTrack(*otherEndNode,track))
 540     {
 541       ((const_cast<Edge *>(this))->_end)->decrRef();
 542       ((const_cast<Edge *>(this))->_end)=otherEndNode;
 543       otherEndNode->incrRef();
 544       return true;
 545     }
 546   return false;
 547 }
 548
 549 /*!
 550  * Precondition : 'start' and 'end' are lying on the same curve than 'this'.
 551  * Add in vec the sub edge lying on this.
 552  * If 'start' is equal (by pointer) to '_end' and 'end' is equal to '_end' too nothing is added.
 553  * If 'start' is equal (by pointer) to '_start' and 'end' is equal to '_start' too nothing is added.
 554  * If 'start' is equal (by pointer) to '_start' and 'end' is equal to '_end' this is added in vec.
 555  */
 556 void Edge::addSubEdgeInVector(Node *start, Node *end, ComposedEdge& vec) const
 557 {
 558   if((start==_start && end==_start) || (start==_end && end==_end))
 559     return ;
 560   if(start==_start && end==_end)
 561     {
 562       incrRef();
 563       vec.pushBack(const_cast<Edge *>(this));
 564       return ;
 565     }
 566   vec.pushBack(buildEdgeLyingOnMe(start,end,true));
 567 }
 568
 569 /*!
 570  * Retrieves a vector 'vectOutput' that is normal to 'this'. 'vectOutput' is normalized.
 571  */
 572 void Edge::getNormalVector(double *vectOutput) const
 573 {
 574   std::copy((const double *)(*_end),(const double *)(*_end)+2,vectOutput);
 575   std::transform(vectOutput,vectOutput+2,(const double *)(*_start),vectOutput,std::minus<double>());
 576   double norm=1./Node::norm(vectOutput);
 577   std::transform(vectOutput,vectOutput+2,vectOutput,bind2nd(std::multiplies<double>(),norm));
 578   double tmp=vectOutput[0];
 579   vectOutput[0]=vectOutput[1];
 580   vectOutput[1]=-tmp;
 581 }
 582
 583 Edge *Edge::BuildEdgeFrom3Points(const double *start, const double *middle, const double *end)
 584 {
 585   Node *b(new Node(start[0],start[1])),*m(new Node(middle[0],middle[1])),*e(new Node(end[0],end[1]));
 586   EdgeLin *e1(new EdgeLin(b,m)),*e2(new EdgeLin(m,e));
 587   SegSegIntersector inters(*e1,*e2); bool colinearity=inters.areColinears(); delete e1; delete e2;
 588   Edge *ret=0;
 589   if(colinearity)
 590     ret=new EdgeLin(b,e);
 591   else
 592     ret=new EdgeArcCircle(b,m,e);
 593   b->decrRef(); m->decrRef(); e->decrRef();
 594   return ret;
 595 }
 596
 597 Edge *Edge::BuildEdgeFrom(Node *start, Node *end)
 598 {
 599   return new EdgeLin(start,end);
 600 }
 601
 602 Edge *Edge::BuildFromXfigLine(std::istream& str)
 603 {
 604   unsigned char type;
 605   str >> type;
 606   if(type=='2')
 607     return new EdgeLin(str);
 608   else if(type=='5')
 609     return new EdgeArcCircle(str);
 610   else
 611     {
 612       std::cerr << "Unknown line found...";
 613       return 0;
 614     }
 615 }
 616
 617 /*!
 618  * \param other The Edge with which we are going to intersect.
 619  * \param commonNode Output. The common nodes found during operation of intersecting.
 620  * \param outVal1 Output filled in case true is returned. It specifies the new or not new edges by which 'this' is replaced after intersecting op.
 621  * \param outVal2 Output filled in case true is returned. It specifies the new or not new edges by which 'other' is replaced after intersecting op.
 622  * return true if the intersection between this.
 623  */
 624 bool Edge::intersectWith(const Edge *other, MergePoints& commonNode,
 625                          ComposedEdge& outVal1, ComposedEdge& outVal2) const
 626 {
 627   bool ret=true;
 628   Bounds *merge=_bounds.nearlyAmIIntersectingWith(other->getBounds());
 629   if(!merge)
 630     return false;
 631   delete merge;
 632   merge=0;
 633   EdgeIntersector *intersector=BuildIntersectorWith(this,other);
 634   ret=Intersect(this,other,intersector,merge,commonNode,outVal1,outVal2);
 635   delete intersector;
 636   return ret;
 637 }
 638
 639 bool Edge::IntersectOverlapped(const Edge *f1, const Edge *f2, EdgeIntersector *intersector, MergePoints& commonNode,
 640                                ComposedEdge& outValForF1, ComposedEdge& outValForF2)
 641 {
 642   bool rev=intersector->haveTheySameDirection();
 643   Node *f2Start=f2->getNode(rev?START:END);
 644   Node *f2End=f2->getNode(rev?END:START);
 645   TypeOfLocInEdge place1, place2;
 646   intersector->getPlacements(f2Start,f2End,place1,place2,commonNode);
 647   int codeForIntersectionCase=CombineCodes(place1,place2);
 648   return SplitOverlappedEdges(f1,f2,f2Start,f2End,rev,codeForIntersectionCase,outValForF1,outValForF2);
 649 }
 650
 651 /*!
 652  * Perform 1D linear interpolation. Warning distrib1 and distrib2 are expected to be in ascending mode.
 653  */
 654 void Edge::Interpolate1DLin(const std::vector<double>& distrib1, const std::vector<double>& distrib2, std::map<int, std::map<int,double> >& result)
 655 {
 656   int nbOfV1=distrib1.size()-1;
 657   int nbOfV2=distrib2.size()-1;
 658   Node *n1=new Node(0.,0.); Node *n3=new Node(0.,0.);
 659   Node *n2=new Node(0.,0.); Node *n4=new Node(0.,0.);
 660   MergePoints commonNode;
 661   for(int i=0;i<nbOfV1;i++)
 662     {
 663       std::vector<double>::const_iterator iter=find_if(distrib2.begin()+1,distrib2.end(),bind2nd(std::greater_equal<double>(),distrib1[i]));
 664       if(iter!=distrib2.end())
 665         {
 666           for(int j=(iter-1)-distrib2.begin();j<nbOfV2;j++)
 667             {
 668               if(distrib2[j]<=distrib1[i+1])
 669                 {
 670                   EdgeLin *e1=new EdgeLin(n1,n2); EdgeLin *e2=new EdgeLin(n3,n4);
 671                   n1->setNewCoords(distrib1[i],0.); n2->setNewCoords(distrib1[i+1],0.);
 672                   n3->setNewCoords(distrib2[j],0.); n4->setNewCoords(distrib2[j+1],0.);
 673                   ComposedEdge *f1=new ComposedEdge;
 674                   ComposedEdge *f2=new ComposedEdge;
 675                   SegSegIntersector inters(*e1,*e2);
 676                   bool b1,b2;
 677                   inters.areOverlappedOrOnlyColinears(0,b1,b2);
 678                   if(IntersectOverlapped(e1,e2,&inters,commonNode,*f1,*f2))
 679                     {
 680                       result[i][j]=f1->getCommonLengthWith(*f2)/e1->getCurveLength();
 681                     }
 682                   ComposedEdge::Delete(f1); ComposedEdge::Delete(f2);
 683                   e1->decrRef(); e2->decrRef();
 684                 }
 685             }
 686         }
 687     }
 688   n1->decrRef(); n2->decrRef(); n3->decrRef(); n4->decrRef();
 689 }
 690
 691 EdgeIntersector *Edge::BuildIntersectorWith(const Edge *e1, const Edge *e2)
 692 {
 693   EdgeIntersector *ret=0;
 694   const EdgeLin *tmp1=0;
 695   const EdgeArcCircle *tmp2=0;
 696   unsigned char type1=e1->getTypeOfFunc();
 697   e1->dynCastFunction(tmp1,tmp2);
 698   unsigned char type2=e2->getTypeOfFunc();
 699   e2->dynCastFunction(tmp1,tmp2);
 700   type1|=type2;
 701   switch(type1)
 702   {
 703     case 1:// Intersection seg/seg
 704       ret=new SegSegIntersector((const EdgeLin &)(*e1),(const EdgeLin &)(*e2));
 705       break;
 706     case 5:// Intersection seg/arc of circle
 707       ret=new ArcCSegIntersector(*tmp2,*tmp1,tmp2==e1);
 708       break;
 709     case 4:// Intersection arc/arc of circle
 710       ret=new ArcCArcCIntersector((const EdgeArcCircle &)(*e1),(const EdgeArcCircle &)(*e2));
 711       break;
 712     default:
 713       //Should never happen
 714       throw Exception("A non managed association of edge has been detected. Go work for intersection computation implementation.");
 715   }
 716   return ret;
 717 }
 718
 719 /*!
 720  * See Node::applySimilarity to see signification of params.
 721  */
 722 void Edge::applySimilarity(double xBary, double yBary, double dimChar)
 723 {
 724   _bounds.applySimilarity(xBary,yBary,dimChar);
 725 }
 726
 727 void Edge::unApplySimilarity(double xBary, double yBary, double dimChar)
 728 {
 729   _bounds.unApplySimilarity(xBary,yBary,dimChar);
 730 }
 731
 732 void Edge::getMiddleOfPointsOriented(const double *p1, const double *p2, double *mid) const
 733 {
 734   return getMiddleOfPoints(p1, p2, mid);
 735 }
 736
 737 bool Edge::Intersect(const Edge *f1, const Edge *f2, EdgeIntersector *intersector, const Bounds *whereToFind, MergePoints& commonNode,
 738                      ComposedEdge& outValForF1, ComposedEdge& outValForF2)
 739 {
 740   bool obviousNoIntersection;
 741   bool areOverlapped;
 742   intersector->areOverlappedOrOnlyColinears(whereToFind,obviousNoIntersection,areOverlapped);
 743   if(areOverlapped)
 744     return IntersectOverlapped(f1,f2,intersector,commonNode,outValForF1,outValForF2);
 745   if(obviousNoIntersection)
 746     return false;
 747   std::vector<Node *> newNodes;
 748   bool order;
 749   if(intersector->intersect(whereToFind,newNodes,order,commonNode))
 750     {
 751       if(newNodes.empty())
 752         throw Exception("Internal error occurred - error in intersector implementation!");// This case should never happen
 753       std::vector<Node *>::iterator iter=newNodes.begin();
 754       std::vector<Node *>::reverse_iterator iterR=newNodes.rbegin();
 755       f1->addSubEdgeInVector(f1->getStartNode(),*iter,outValForF1);
 756       f2->addSubEdgeInVector(f2->getStartNode(),order?*iter:*iterR,outValForF2);
 757       for(std::vector<Node *>::iterator iter2=newNodes.begin();iter2!=newNodes.end();iter2++,iterR++)
 758         {
 759           if((iter2+1)==newNodes.end())
 760             {
 761               f1->addSubEdgeInVector(*iter2,f1->getEndNode(),outValForF1);
 762               (*iter2)->decrRef();
 763               f2->addSubEdgeInVector(order?*iter2:*iterR,f2->getEndNode(),outValForF2);
 764             }
 765           else
 766             {
 767               f1->addSubEdgeInVector(*iter2,*(iter2+1),outValForF1);
 768               (*iter2)->decrRef();
 769               f2->addSubEdgeInVector(order?*iter2:*iterR,order?*(iter2+1):*(iterR+1),outValForF2);
 770             }
 771         }
 772       return true;
 773     }
 774   else//no intersection inside whereToFind
 775     return false;
 776 }
 777
 778 int Edge::CombineCodes(TypeOfLocInEdge code1, TypeOfLocInEdge code2)
 779 {
 780   int ret=(int)code1;
 781   ret*=OFFSET_FOR_TYPEOFLOCINEDGE;
 782   ret+=(int)code2;
 783   return ret;
 784 }
 785
 786 /*!
 787  * This method splits e1 and e2 into pieces as much sharable as possible. The precondition to the call of this method
 788  * is that e1 and e2 have been declared as overlapped by corresponding intersector built from e1 and e2 type.
 789  *
 790  * @param nS start node of e2 with the SAME DIRECTION as e1. The pointer nS should be equal to start node of e2 or to its end node.
 791  * @param nE end node of e2 with the SAME DIRECTION as e1. The pointer nE should be equal to start node of e2 or to its end node.
 792  * @param direction is param that specifies if e2 and e1 have same directions (true) or opposed (false).
 793  * @param code is the code returned by method Edge::combineCodes.
 794  */
 795 bool Edge::SplitOverlappedEdges(const Edge *e1, const Edge *e2, Node *nS, Node *nE, bool direction, int code,
 796                                 ComposedEdge& outVal1, ComposedEdge& outVal2)
 797 {
 798   Edge *tmp;
 799   switch(code)
 800   {
 801     case OUT_BEFORE*OFFSET_FOR_TYPEOFLOCINEDGE+START:      // OUT_BEFORE - START
 802     case OUT_BEFORE*OFFSET_FOR_TYPEOFLOCINEDGE+OUT_BEFORE: // OUT_BEFORE - OUT_BEFORE
 803     case OUT_AFTER*OFFSET_FOR_TYPEOFLOCINEDGE+OUT_AFTER:   // OUT_AFTER - OUT_AFTER
 804     case END*OFFSET_FOR_TYPEOFLOCINEDGE+OUT_AFTER:         // END - OUT_AFTER
 805     case END*OFFSET_FOR_TYPEOFLOCINEDGE+START:             // END - START
 806     return false;
 807     case INSIDE*OFFSET_FOR_TYPEOFLOCINEDGE+OUT_AFTER:      // INSIDE - OUT_AFTER
 808     outVal1.pushBack(e1->buildEdgeLyingOnMe(e1->getStartNode(),nS,true));
 809     tmp=e1->buildEdgeLyingOnMe(nS,e1->getEndNode()); tmp->incrRef();
 810     outVal1.pushBack(tmp);
 811     outVal2.resize(2);
 812     outVal2.setValueAt(direction?0:1,tmp,direction); tmp->declareOn();
 813     outVal2.setValueAt(direction?1:0,e1->buildEdgeLyingOnMe(e1->getEndNode(),nE,direction));
 814     return true;
 815     case INSIDE*OFFSET_FOR_TYPEOFLOCINEDGE+INSIDE:         // INSIDE - INSIDE
 816     {
 817       if(!e2->isIn(e2->getCharactValue(*(e1->getStartNode()))))
 818         {
 819           e2->incrRef(); e2->incrRef();
 820           outVal1.resize(3);
 821           outVal1.setValueAt(0,e1->buildEdgeLyingOnMe(e1->getStartNode(),nS));
 822           outVal1.setValueAt(1,const_cast<Edge*>(e2),direction);
 823           outVal1.setValueAt(2,e1->buildEdgeLyingOnMe(nE,e1->getEndNode()));
 824           outVal2.pushBack(const_cast<Edge*>(e2)); e2->declareOn();
 825           return true;
 826         }
 827       else
 828         {
 829           outVal1.resize(3);
 830           outVal2.resize(3);
 831           tmp=e1->buildEdgeLyingOnMe(e1->getStartNode(),nE); tmp->incrRef(); tmp->declareOn();
 832           outVal1.setValueAt(0,tmp,true); outVal2.setValueAt(direction?2:0,tmp,direction);
 833           outVal1.setValueAt(1,e1->buildEdgeLyingOnMe(nE,nS));
 834           tmp=e1->buildEdgeLyingOnMe(nS,e1->getEndNode()); tmp->incrRef(); tmp->declareOn();
 835           outVal1.setValueAt(2,tmp,true); outVal2.setValueAt(direction?0:2,tmp,direction);
 836           tmp=e1->buildEdgeLyingOnMe(e1->getEndNode(),e1->getStartNode());
 837           outVal2.setValueAt(1,tmp,direction);
 838           return true;
 839         }
 840     }
 841     case OUT_BEFORE*OFFSET_FOR_TYPEOFLOCINEDGE+INSIDE:     // OUT_BEFORE - INSIDE
 842     {
 843       tmp=e1->buildEdgeLyingOnMe(e1->getStartNode(),nE); tmp->incrRef();
 844       outVal1.pushBack(tmp);
 845       outVal1.pushBack(e1->buildEdgeLyingOnMe(nE,e1->getEndNode()));
 846       outVal2.resize(2);
 847       outVal2.setValueAt(direction?0:1,e1->buildEdgeLyingOnMe(nS,e1->getStartNode(),direction));
 848       outVal2.setValueAt(direction?1:0,tmp,direction); tmp->declareOn();
 849       return true;
 850     }
 851     case OUT_BEFORE*OFFSET_FOR_TYPEOFLOCINEDGE+OUT_AFTER:  // OUT_BEFORE - OUT_AFTER
 852     {
 853       e1->incrRef(); e1->incrRef();
 854       outVal1.pushBack(const_cast<Edge*>(e1));
 855       outVal2.resize(3);
 856       outVal2.setValueAt(direction?0:2,e1->buildEdgeLyingOnMe(nS,e1->getStartNode(),direction));
 857       outVal2.setValueAt(1,const_cast<Edge*>(e1),direction); e1->declareOn();
 858       outVal2.setValueAt(direction?2:0,e1->buildEdgeLyingOnMe(e1->getEndNode(),nE,direction));
 859       return true;
 860     }
 861     case START*OFFSET_FOR_TYPEOFLOCINEDGE+END:             // START - END
 862     {
 863       e1->incrRef(); e1->incrRef();
 864       outVal1.pushBack(const_cast<Edge*>(e1));
 865       outVal2.pushBack(const_cast<Edge*>(e1),direction); e1->declareOn();
 866       return true;
 867     }
 868     case START*OFFSET_FOR_TYPEOFLOCINEDGE+OUT_AFTER:       // START - OUT_AFTER
 869     {
 870       e1->incrRef(); e1->incrRef();
 871       outVal1.pushBack(const_cast<Edge*>(e1));
 872       outVal2.resize(2);
 873       outVal2.setValueAt(direction?0:1,const_cast<Edge*>(e1),direction); e1->declareOn();
 874       outVal2.setValueAt(direction?1:0,e1->buildEdgeLyingOnMe(e1->getEndNode(),nE,direction));
 875       return true;
 876     }
 877     case INSIDE*OFFSET_FOR_TYPEOFLOCINEDGE+END:            // INSIDE - END
 878     {
 879       e2->incrRef(); e2->incrRef();
 880       outVal1.pushBack(e1->buildEdgeLyingOnMe(e1->getStartNode(),nS,true));
 881       outVal1.pushBack(const_cast<Edge*>(e2),direction);
 882       outVal2.pushBack(const_cast<Edge*>(e2)); e2->declareOn();
 883       return true;
 884     }
 885     case OUT_BEFORE*OFFSET_FOR_TYPEOFLOCINEDGE+END:        // OUT_BEFORE - END
 886     {
 887       e1->incrRef(); e1->incrRef();
 888       outVal1.pushBack(const_cast<Edge*>(e1));
 889       outVal2.resize(2);
 890       outVal2.setValueAt(direction?0:1,e1->buildEdgeLyingOnMe(nS,e1->getStartNode(),direction));
 891       outVal2.setValueAt(direction?1:0,const_cast<Edge*>(e1),direction); e1->declareOn();
 892       return true;
 893     }
 894     case START*OFFSET_FOR_TYPEOFLOCINEDGE+INSIDE:          // START - INSIDE
 895     {
 896       e2->incrRef(); e2->incrRef();
 897       outVal1.pushBack(const_cast<Edge*>(e2),direction);
 898       outVal1.pushBack(e1->buildEdgeLyingOnMe(nE,e1->getEndNode()));
 899       outVal2.pushBack(const_cast<Edge*>(e2)); e2->declareOn();
 900       return true;
 901     }
 902     case INSIDE*OFFSET_FOR_TYPEOFLOCINEDGE+START:          // INSIDE - START
 903     {
 904       outVal1.resize(2);
 905       outVal2.resize(2);
 906       tmp=e1->buildEdgeLyingOnMe(nS,e1->getEndNode()); tmp->incrRef(); tmp->declareOn();
 907       outVal1.setValueAt(0,e1->buildEdgeLyingOnMe(e1->getStartNode(),nS));
 908       outVal1.setValueAt(1,tmp);
 909       outVal2.setValueAt(direction?0:1,tmp,direction);
 910       outVal2.setValueAt(direction?1:0,e1->buildEdgeLyingOnMe(e1->getEndNode(),nE,direction));
 911       return true;
 912     }
 913     case END*OFFSET_FOR_TYPEOFLOCINEDGE+INSIDE:            // END - INSIDE
 914     {
 915       outVal1.resize(2);
 916       outVal2.resize(2);
 917       tmp=e1->buildEdgeLyingOnMe(e1->getStartNode(),nE); tmp->incrRef(); tmp->declareOn();
 918       outVal1.setValueAt(0,tmp);
 919       outVal1.setValueAt(1,e1->buildEdgeLyingOnMe(nE,e1->getEndNode()));
 920       outVal2.setValueAt(direction?0:1,e1->buildEdgeLyingOnMe(e1->getEndNode(),e1->getStartNode(),direction));
 921       outVal2.setValueAt(direction?1:0,tmp,direction);
 922       return true;
 923     }
 924     default:
 925       throw Exception("Unexpected situation of overlapping edges : internal error occurs ! ");
 926   }
 927 }
 928
 929 void Edge::dumpToCout(const std::map<INTERP_KERNEL::Node *,int>& mapp, int index) const
 930 {
 931   auto sI(mapp.find(getStartNode())), eI(mapp.find(getEndNode()));
 932   int start = (sI == mapp.end() ? -1 : sI->second), end = (eI == mapp.end() ? -1 : eI->second);
 933   std::string locs;
 934   switch (getLoc())
 935   {
 936     case FULL_IN_1: locs="FULL_IN_1"; break;
 937     case FULL_ON_1: locs="FULL_ON_1"; break;
 938     case FULL_OUT_1: locs="FULL_OUT_1"; break;
 939     case FULL_UNKNOWN: locs="FULL_UNKNOWN"; break;
 940     default: locs="oh my God! This is so wrong.";
 941   }
 942   std::cout << "Edge [" << index << "] : ("<<  std::hex << this << std::dec << ") -> (" << start << ", " << end << ")\t" << locs << std::endl;
 943 }
 944
 945 bool Edge::isEqual(const Edge& other) const
 946 {
 947   return _start->isEqual(*other._start) && _end->isEqual(*other._end);
 948 }
 949
 950 /**
 951  * This method takes in input nodes in \a subNodes (using \a coo)
 952  *
 953  * \param [in,out] subNodes to be sorted
 954  * \return true if a reordering was necessary false if not.
 955  */
 956 bool Edge::sortSubNodesAbs(const double *coo, std::vector<int>& subNodes)
 957 {
 958   Bounds b;
 959   b.prepareForAggregation();
 960   b.aggregate(getBounds());
 961   double xBary,yBary;
 962   double dimChar(b.getCaracteristicDim());
 963   b.getBarycenter(xBary,yBary);
 964   applySimilarity(xBary,yBary,dimChar);
 965   _start->applySimilarity(xBary,yBary,dimChar);
 966   _end->applySimilarity(xBary,yBary,dimChar);
 967   //
 968   std::size_t sz(subNodes.size()),i(0);
 969   std::vector< std::pair<double,Node *> > an2(sz);
 970   std::map<Node *, int> m;
 971   for(std::vector<int>::const_iterator it=subNodes.begin();it!=subNodes.end();it++,i++)
 972     {
 973       Node *n(new Node(coo[2*(*it)],coo[2*(*it)+1]));
 974       n->applySimilarity(xBary,yBary,dimChar);
 975       m[n]=*it;
 976       an2[i]=std::pair<double,Node *>(getCharactValueBtw0And1(*n),n);
 977     }
 978   std::sort(an2.begin(),an2.end());
 979   //
 980   bool ret(false);
 981   for(i=0;i<sz;i++)
 982     {
 983       int id(m[an2[i].second]);
 984       if(id!=subNodes[i])
 985         { subNodes[i]=id; ret=true; }
 986     }
 987   //
 988   for(std::map<INTERP_KERNEL::Node *,int>::const_iterator it2=m.begin();it2!=m.end();it2++)
 989     (*it2).first->decrRef();
 990   return ret;
 991 }
 992
 993 /**
 994  * Sort nodes so that they all lie consecutively on the edge that has been cut.
 995  */
 996 void Edge::sortIdsAbs(const std::vector<INTERP_KERNEL::Node *>& addNodes, const std::map<INTERP_KERNEL::Node *, int>& mapp1,
 997                       const std::map<INTERP_KERNEL::Node *, int>& mapp2, std::vector<int>& edgesThis)
 998 {
 999   int startId=(*mapp1.find(_start)).second;
1000   int endId=(*mapp1.find(_end)).second;
1001   if (! addNodes.size()) // quick way out, no new node to add.
1002     {
1003       edgesThis.push_back(startId);
1004       edgesThis.push_back(endId);
1005       return;
1006     }
1007
1008   Bounds b;
1009   b.prepareForAggregation();
1010   b.aggregate(getBounds());
1011   double xBary,yBary;
1012   double dimChar=b.getCaracteristicDim();
1013   b.getBarycenter(xBary,yBary);
1014   for(std::vector<Node *>::const_iterator iter=addNodes.begin();iter!=addNodes.end();iter++)
1015     (*iter)->applySimilarity(xBary,yBary,dimChar);
1016   applySimilarity(xBary,yBary,dimChar);
1017   _start->applySimilarity(xBary,yBary,dimChar);
1018   _end->applySimilarity(xBary,yBary,dimChar);
1019   std::size_t sz=addNodes.size();
1020   std::vector< std::pair<double,Node *> > an2(sz);
1021   for(std::size_t i=0;i<sz;i++)
1022     an2[i]=std::pair<double,Node *>(getCharactValueBtw0And1(*addNodes[i]),addNodes[i]);
1023   std::sort(an2.begin(),an2.end());
1024   std::vector<int> tmpp;
1025   for(std::vector< std::pair<double,Node *> >::const_iterator it=an2.begin();it!=an2.end();it++)
1026     {
1027       int idd=(*mapp2.find((*it).second)).second;
1028       if((*it).first<QuadraticPlanarPrecision::getPrecision())
1029         {
1030           startId=idd;
1031           continue;
1032         }
1033       if((*it).first>1-QuadraticPlanarPrecision::getPrecision())
1034         {
1035           endId=idd;
1036           continue;
1037         }
1038       tmpp.push_back(idd);
1039     }
1040   std::vector<int> tmpp2(tmpp.size()+2);
1041   tmpp2[0]=startId;
1042   std::copy(tmpp.begin(),tmpp.end(),tmpp2.begin()+1);
1043   tmpp2[tmpp.size()+1]=endId;
1044   std::vector<int>::iterator itt=std::unique(tmpp2.begin(),tmpp2.end());
1045   tmpp2.resize(std::distance(tmpp2.begin(),itt));
1046   int nbOfEdges=tmpp2.size()-1;
1047   for(int i=0;i<nbOfEdges;i++)
1048     {
1049       edgesThis.push_back(tmpp2[i]);
1050       edgesThis.push_back(tmpp2[i+1]);
1051     }
1052 }
1053
1054 void Edge::fillGlobalInfoAbs(bool direction, const std::map<INTERP_KERNEL::Node *,int>& mapThis, const std::map<INTERP_KERNEL::Node *,int>& mapOther, int offset1, int offset2, double fact, double baryX, double baryY,
1055                                 std::vector<int>& edgesThis, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int> mapAddCoo) const
1056 {
1057   int tmp[2];
1058   _start->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp);
1059   _end->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,tmp+1);
1060   if(direction)
1061     {
1062       edgesThis.push_back(tmp[0]);
1063       edgesThis.push_back(tmp[1]);
1064     }
1065   else
1066     {
1067       edgesThis.push_back(tmp[1]);
1068       edgesThis.push_back(tmp[0]);
1069     }
1070 }
1071
1072 void Edge::fillGlobalInfoAbs2(const std::map<INTERP_KERNEL::Node *,int>& mapThis, const std::map<INTERP_KERNEL::Node *,int>& mapOther, int offset1, int offset2, double fact, double baryX, double baryY,
1073                               unsigned skipStartOrEnd,
1074                               std::vector<int>& edgesOther, std::vector<double>& addCoo, std::map<INTERP_KERNEL::Node *,int>& mapAddCoo) const
1075 {
1076   if (skipStartOrEnd != 1) // see meaning in splitAbs()
1077     _start->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther);
1078   if (skipStartOrEnd != 2)
1079   _end->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,baryX,baryY,addCoo,mapAddCoo,edgesOther);
1080 }