src/INTERP_KERNEL/Geometric2D/InterpKernelGeo2DQuadraticPolygon.cxx

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