Back to index

salome-med  6.5.0
InterpKernelGeo2DQuadraticPolygon.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "InterpKernelGeo2DQuadraticPolygon.hxx"
00021 #include "InterpKernelGeo2DElementaryEdge.hxx"
00022 #include "InterpKernelGeo2DEdgeArcCircle.hxx"
00023 #include "InterpKernelGeo2DAbstractEdge.hxx"
00024 #include "InterpKernelGeo2DEdgeLin.hxx"
00025 #include "InterpKernelGeo2DBounds.hxx"
00026 #include "InterpKernelGeo2DEdge.txx"
00027 
00028 #include "NormalizedUnstructuredMesh.hxx"
00029 
00030 #include <fstream>
00031 #include <iomanip>
00032 #include <cstring>
00033 #include <limits>
00034 
00035 using namespace INTERP_KERNEL;
00036 
00037 namespace INTERP_KERNEL
00038 {
00039   const unsigned MAX_SIZE_OF_LINE_XFIG_FILE=1024;
00040 }
00041 
00042 QuadraticPolygon::QuadraticPolygon(const char *file)
00043 {
00044   char currentLine[MAX_SIZE_OF_LINE_XFIG_FILE];
00045   std::ifstream stream(file);
00046   stream.exceptions(std::ios_base::eofbit);
00047   try
00048     {
00049       do
00050         stream.getline(currentLine,MAX_SIZE_OF_LINE_XFIG_FILE);
00051       while(strcmp(currentLine,"1200 2")!=0);
00052       do
00053         {
00054           Edge *newEdge=Edge::BuildFromXfigLine(stream);
00055           if(!empty())
00056             newEdge->changeStartNodeWith(back()->getEndNode());
00057           pushBack(newEdge);
00058         }
00059       while(1);
00060     }
00061   catch(std::ifstream::failure&)
00062     {
00063     }
00064   front()->changeStartNodeWith(back()->getEndNode());
00065 }
00066 
00067 QuadraticPolygon::~QuadraticPolygon()
00068 {
00069 }
00070 
00071 QuadraticPolygon *QuadraticPolygon::BuildLinearPolygon(std::vector<Node *>& nodes)
00072 {
00073   QuadraticPolygon *ret=new QuadraticPolygon;
00074   std::size_t size=nodes.size();
00075   for(std::size_t i=0;i<size;i++)
00076     {
00077       ret->pushBack(new EdgeLin(nodes[i],nodes[(i+1)%size]));
00078       nodes[i]->decrRef();
00079     }
00080   return ret;
00081 }
00082 
00083 QuadraticPolygon *QuadraticPolygon::BuildArcCirclePolygon(std::vector<Node *>& nodes)
00084 {
00085   QuadraticPolygon *ret=new QuadraticPolygon;
00086   std::size_t size=nodes.size();
00087   for(std::size_t i=0;i<size/2;i++)
00088     {
00089       EdgeLin *e1,*e2;
00090       e1=new EdgeLin(nodes[i],nodes[i+size/2]);
00091       e2=new EdgeLin(nodes[i+size/2],nodes[(i+1)%(size/2)]);
00092       SegSegIntersector inters(*e1,*e2);
00093       bool colinearity=inters.areColinears();
00094       delete e1; delete e2;
00095       if(colinearity)
00096         ret->pushBack(new EdgeLin(nodes[i],nodes[(i+1)%(size/2)]));
00097       else
00098         ret->pushBack(new EdgeArcCircle(nodes[i],nodes[i+size/2],nodes[(i+1)%(size/2)]));
00099       nodes[i]->decrRef(); nodes[i+size/2]->decrRef();
00100     }
00101   return ret;
00102 }
00103 
00104 void QuadraticPolygon::BuildDbgFile(const std::vector<Node *>& nodes, const char *fileName)
00105 {
00106   std::ofstream file(fileName);
00107   file << std::setprecision(16);
00108   file << "  double coords[]=" << std::endl << "    { ";
00109   for(std::vector<Node *>::const_iterator iter=nodes.begin();iter!=nodes.end();iter++)
00110     {
00111       if(iter!=nodes.begin())
00112         file << "," << std::endl << "      ";
00113       file << (*(*iter))[0] << ", " << (*(*iter))[1];
00114     }
00115   file << "};" << std::endl;
00116 }
00117 
00118 void QuadraticPolygon::closeMe() const
00119 {
00120   if(!front()->changeStartNodeWith(back()->getEndNode()))
00121     throw(Exception("big error: not closed polygon..."));
00122 }
00123 
00124 void QuadraticPolygon::circularPermute()
00125 {
00126   if(_sub_edges.size()>1)
00127     {
00128       ElementaryEdge *first=_sub_edges.front();
00129       _sub_edges.pop_front();
00130       _sub_edges.push_back(first);
00131     }
00132 }
00133 
00134 bool QuadraticPolygon::isButterflyAbs()
00135 {
00136   INTERP_KERNEL::Bounds b;
00137   double xBary,yBary;
00138   b.prepareForAggregation();
00139   fillBounds(b); 
00140   double dimChar=b.getCaracteristicDim();
00141   b.getBarycenter(xBary,yBary);
00142   applyGlobalSimilarity(xBary,yBary,dimChar);
00143   //
00144   return isButterfly();
00145 }
00146 
00147 bool QuadraticPolygon::isButterfly() const
00148 {
00149   for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++)
00150     {
00151       Edge *e1=(*it)->getPtr();
00152       std::list<ElementaryEdge *>::const_iterator it2=it;
00153       it2++;
00154       for(;it2!=_sub_edges.end();it2++)
00155         {
00156           MergePoints commonNode;
00157           ComposedEdge *outVal1=new ComposedEdge;
00158           ComposedEdge *outVal2=new ComposedEdge;
00159           Edge *e2=(*it2)->getPtr();
00160           if(e1->intersectWith(e2,commonNode,*outVal1,*outVal2))
00161             {
00162               Delete(outVal1);
00163               Delete(outVal2);
00164               return true;
00165             }
00166           Delete(outVal1);
00167           Delete(outVal2);
00168         }
00169     }
00170   return false;
00171 }
00172 
00173 void QuadraticPolygon::dumpInXfigFileWithOther(const ComposedEdge& other, const char *fileName) const
00174 {
00175   std::ofstream file(fileName);
00176   const int resolution=1200;
00177   Bounds box;
00178   box.prepareForAggregation();
00179   fillBounds(box);
00180   other.fillBounds(box);
00181   dumpInXfigFile(file,resolution,box);
00182   other.ComposedEdge::dumpInXfigFile(file,resolution,box);
00183 }
00184 
00185 void QuadraticPolygon::dumpInXfigFile(const char *fileName) const
00186 {
00187   std::ofstream file(fileName);
00188   const int resolution=1200;
00189   Bounds box;
00190   box.prepareForAggregation();
00191   fillBounds(box);
00192   dumpInXfigFile(file,resolution,box);
00193 }
00194 
00195 void QuadraticPolygon::dumpInXfigFile(std::ostream& stream, int resolution, const Bounds& box) const
00196 {
00197   stream << "#FIG 3.2  Produced by xfig version 3.2.5-alpha5" << std::endl;
00198   stream << "Landscape" << std::endl;
00199   stream << "Center" << std::endl;
00200   stream << "Metric" << std::endl;
00201   stream << "Letter" << std::endl;
00202   stream << "100.00" << std::endl;
00203   stream << "Single" << std::endl;
00204   stream << "-2" << std::endl;
00205   stream << resolution << " 2" << std::endl;
00206   ComposedEdge::dumpInXfigFile(stream,resolution,box);
00207 }
00208 
00212 double QuadraticPolygon::intersectWithAbs(QuadraticPolygon& other)
00213 {
00214   double ret=0.,xBaryBB,yBaryBB;
00215   double fact=normalize(&other,xBaryBB,yBaryBB);
00216   std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
00217   for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
00218     {
00219       ret+=fabs((*iter)->getArea());
00220       delete *iter;
00221     }
00222   return ret*fact*fact;
00223 }
00224 
00236 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,
00237                                 std::vector<int>& edgesThis, int cellIdThis, std::vector< std::vector<int> >& edgesInOtherColinearWithThis, std::vector< std::vector<int> >& subDivOther, std::vector<double>& addCoo)
00238 {
00239   double xBaryBB, yBaryBB;
00240   double fact=normalizeExt(&other, xBaryBB, yBaryBB);
00241   //
00242   IteratorOnComposedEdge it1(this),it3(&other);
00243   MergePoints merge;
00244   ComposedEdge *c1=new ComposedEdge;
00245   ComposedEdge *c2=new ComposedEdge;
00246   int i=0;
00247   std::map<INTERP_KERNEL::Node *,int> mapAddCoo;
00248   for(it3.first();!it3.finished();it3.next(),i++)//iteration over 'other' _sub_edges
00249     {
00250       QuadraticPolygon otherTmp;
00251       ElementaryEdge* curE3=it3.current();
00252       otherTmp.pushBack(new ElementaryEdge(curE3->getPtr(),curE3->getDirection())); curE3->getPtr()->incrRef();
00253       IteratorOnComposedEdge it2(&otherTmp);
00254       for(it2.first();!it2.finished();it2.next())//iteration on subedges of 'other->_sub_edge'
00255         {
00256           ElementaryEdge* curE2=it2.current();
00257           if(!curE2->isThereStartPoint())
00258             it1.first();
00259           else
00260             it1=curE2->getIterator();
00261           for(;!it1.finished();)//iteration over 'this' _sub_edges
00262             {
00263               ElementaryEdge* curE1=it1.current();
00264               merge.clear();
00265               if(curE1->getPtr()->intersectWith(curE2->getPtr(),merge,*c1,*c2))
00266                 {
00267                   if(!curE1->getDirection()) c1->reverse();
00268                   if(!curE2->getDirection()) c2->reverse();
00269                   UpdateNeighbours(merge,it1,it2,c1,c2);
00270                   //Substitution of simple edge by sub-edges.
00271                   delete curE1; // <-- destroying simple edge coming from pol1
00272                   delete curE2; // <-- destroying simple edge coming from pol2
00273                   it1.insertElemEdges(c1,true);// <-- 2nd param is true to go next.
00274                   it2.insertElemEdges(c2,false);// <-- 2nd param is false to avoid to go next.
00275                   curE2=it2.current();
00276                   //
00277                   it1.assignMySelfToAllElems(c2);//To avoid that others
00278                   SoftDelete(c1);
00279                   SoftDelete(c2);
00280                   c1=new ComposedEdge;
00281                   c2=new ComposedEdge;
00282                 }
00283               else
00284                 {
00285                   UpdateNeighbours(merge,it1,it2,curE1,curE2);
00286                   it1.next();
00287                 }
00288             }
00289         }
00290       if(otherTmp.presenceOfOn())
00291         edgesInOtherColinearWithThis[otherEdgeIds[i]].push_back(cellIdThis);
00292       if(otherTmp._sub_edges.size()>1)
00293         {
00294           for(std::list<ElementaryEdge *>::const_iterator it=otherTmp._sub_edges.begin();it!=otherTmp._sub_edges.end();it++)
00295             (*it)->fillGlobalInfoAbs2(mapThis,mapOther,offset1,offset2,fact,xBaryBB,yBaryBB,subDivOther[otherEdgeIds[i]],addCoo,mapAddCoo);
00296         }
00297     }
00298   Delete(c1);
00299   Delete(c2);
00300   //
00301   for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++)
00302     (*it)->fillGlobalInfoAbs(mapThis,mapOther,offset1,offset2,fact,xBaryBB,yBaryBB,edgesThis,addCoo,mapAddCoo);
00303   //
00304 }
00305 
00310 void QuadraticPolygon::buildFromCrudeDataArray(const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad, const int *nodalBg, const double *coords,
00311                                                const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges)
00312 {
00313   std::size_t nbOfSeg=std::distance(descBg,descEnd);
00314   for(std::size_t i=0;i<nbOfSeg;i++)
00315     {
00316       appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
00317     }
00318 }
00319 
00320 void QuadraticPolygon::appendEdgeFromCrudeDataArray(std::size_t edgePos, const std::map<int,INTERP_KERNEL::Node *>& mapp, bool isQuad, const int *nodalBg, const double *coords,
00321                                                     const int *descBg, const int *descEnd, const std::vector<std::vector<int> >& intersectEdges)
00322 {
00323   if(!isQuad)
00324     {
00325       bool direct=descBg[edgePos]>0;
00326       int edgeId=abs(descBg[edgePos])-1;
00327       const std::vector<int>& subEdge=intersectEdges[edgeId];
00328       std::size_t nbOfSubEdges=subEdge.size()/2;
00329       for(std::size_t j=0;j<nbOfSubEdges;j++)
00330         appendSubEdgeFromCrudeDataArray(0,j,direct,edgeId,subEdge,mapp);
00331     }
00332   else
00333     {
00334       std::size_t nbOfSeg=std::distance(descBg,descEnd);
00335       const double *st=coords+2*(nodalBg[edgePos]); 
00336       INTERP_KERNEL::Node *st0=new INTERP_KERNEL::Node(st[0],st[1]);
00337       const double *endd=coords+2*(nodalBg[(edgePos+1)%nbOfSeg]);
00338       INTERP_KERNEL::Node *endd0=new INTERP_KERNEL::Node(endd[0],endd[1]);
00339       const double *middle=coords+2*(nodalBg[edgePos+nbOfSeg]);
00340       INTERP_KERNEL::Node *middle0=new INTERP_KERNEL::Node(middle[0],middle[1]);
00341       EdgeLin *e1,*e2;
00342       e1=new EdgeLin(st0,middle0);
00343       e2=new EdgeLin(middle0,endd0);
00344       SegSegIntersector inters(*e1,*e2);
00345       bool colinearity=inters.areColinears();
00346       delete e1; delete e2;
00347       //
00348       bool direct=descBg[edgePos]>0;
00349       int edgeId=abs(descBg[edgePos])-1;
00350       const std::vector<int>& subEdge=intersectEdges[edgeId];
00351       std::size_t nbOfSubEdges=subEdge.size()/2;
00352       if(colinearity)
00353         {   
00354           for(std::size_t j=0;j<nbOfSubEdges;j++)
00355             appendSubEdgeFromCrudeDataArray(0,j,direct,edgeId,subEdge,mapp);
00356         }
00357       else
00358         {
00359           Edge *e=new EdgeArcCircle(st0,middle0,endd0,true);
00360           for(std::size_t j=0;j<nbOfSubEdges;j++)
00361             appendSubEdgeFromCrudeDataArray(e,j,direct,edgeId,subEdge,mapp);
00362           e->decrRef();
00363         }
00364       st0->decrRef(); endd0->decrRef(); middle0->decrRef();
00365     }
00366 }
00367 
00368 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)
00369 {
00370   std::size_t nbOfSubEdges=subEdge.size()/2;
00371   if(!baseEdge)
00372     {//it is not a quadratic subedge
00373       Node *start=(*mapp.find(direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1])).second;
00374       Node *end=(*mapp.find(direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2])).second;
00375       ElementaryEdge *e=ElementaryEdge::BuildEdgeFromCrudeDataArray(true,start,end);
00376       pushBack(e);
00377     }
00378   else
00379     {//it is a quadratic subedge
00380       Node *start=(*mapp.find(direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1])).second;
00381       Node *end=(*mapp.find(direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2])).second;
00382       Edge *ee=baseEdge->buildEdgeLyingOnMe(start,end);
00383       ElementaryEdge *eee=new ElementaryEdge(ee,true);
00384       pushBack(eee);
00385     }
00386 }
00387 
00392 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,
00393                                                 const INTERP_KERNEL::QuadraticPolygon& pol1, const int *descBg1, const int *descEnd1, const std::vector<std::vector<int> >& intersectEdges1,
00394                                                 const std::vector< std::vector<int> >& colinear1)
00395 {
00396   std::size_t nbOfSeg=std::distance(descBg,descEnd);
00397   for(std::size_t i=0;i<nbOfSeg;i++)//loop over all edges of pol2
00398     {
00399       bool direct=descBg[i]>0;
00400       int edgeId=abs(descBg[i])-1;//current edge id of pol2
00401       bool directos=colinear1[edgeId].empty();
00402       std::vector<std::pair<int,std::pair<bool,int> > > idIns1;
00403       int offset1=0;
00404       if(!directos)
00405         {// if the current edge of pol2 has one or more colinear edges part into pol1
00406           const std::vector<int>& c=colinear1[edgeId];
00407           std::size_t nbOfEdgesIn1=std::distance(descBg1,descEnd1);
00408           for(std::size_t j=0;j<nbOfEdgesIn1;j++)
00409             {
00410               int edgeId1=abs(descBg1[j])-1;
00411               if(std::find(c.begin(),c.end(),edgeId1)!=c.end())
00412                 {
00413                   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
00414                   //std::pair<edgeId1); direct1=descBg1[j]>0;
00415                 }
00416               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
00417             }
00418           directos=idIns1.empty();
00419         }
00420       if(directos)
00421         {//no subpart of edge 'edgeId' of pol2 is in pol1 so let's operate the same thing that QuadraticPolygon::buildFromCrudeDataArray method
00422           appendEdgeFromCrudeDataArray(i,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
00423         }
00424       else
00425         {//there is subpart of edge 'edgeId' of pol2 inside pol1
00426           const std::vector<int>& subEdge=intersectEdges[edgeId];
00427           std::size_t nbOfSubEdges=subEdge.size()/2;
00428           for(std::size_t j=0;j<nbOfSubEdges;j++)
00429             {
00430               int idBg=direct?subEdge[2*j]:subEdge[2*nbOfSubEdges-2*j-1];
00431               int idEnd=direct?subEdge[2*j+1]:subEdge[2*nbOfSubEdges-2*j-2];
00432               bool direction11,found=false;
00433               bool direct1;//store if needed the direction in 1
00434               int offset2;
00435               std::size_t nbOfSubEdges1;
00436               for(std::vector<std::pair<int,std::pair<bool,int> > >::const_iterator it=idIns1.begin();it!=idIns1.end();it++)
00437                 {
00438                   int idIn1=(*it).first;//store if needed the cell id in 1
00439                   direct1=(*it).second.first;
00440                   offset1=(*it).second.second;
00441                   const std::vector<int>& subEdge1PossiblyAlreadyIn1=intersectEdges1[idIn1];
00442                   nbOfSubEdges1=subEdge1PossiblyAlreadyIn1.size()/2;
00443                   offset2=0;
00444                   for(std::size_t k=0;k<nbOfSubEdges1 && !found;k++)
00445                     {//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
00446                       if(subEdge1PossiblyAlreadyIn1[2*k]==idBg && subEdge1PossiblyAlreadyIn1[2*k+1]==idEnd)
00447                         { direction11=true; found=true; }
00448                       else if(subEdge1PossiblyAlreadyIn1[2*k]==idEnd && subEdge1PossiblyAlreadyIn1[2*k+1]==idBg)
00449                         { direction11=false; found=true; }
00450                       else
00451                         offset2++;
00452                     }
00453                 }
00454               if(!found)
00455                 {//the current subedge of edge 'edgeId' of pol2 is not a part of the colinear edge 'idIn1' of pol1 -> build new Edge instance
00456                   //appendEdgeFromCrudeDataArray(j,mapp,isQuad,nodalBg,coords,descBg,descEnd,intersectEdges);
00457                   Node *start=(*mapp.find(idBg)).second;
00458                   Node *end=(*mapp.find(idEnd)).second;
00459                   ElementaryEdge *e=ElementaryEdge::BuildEdgeFromCrudeDataArray(true,start,end);
00460                   pushBack(e);
00461                 }
00462               else
00463                 {//the current subedge of edge 'edgeId' of pol2 is part of the colinear edge 'idIn1' of pol1 -> reuse Edge instance of pol1
00464                   ElementaryEdge *e=pol1[offset1+(direct1?offset2:nbOfSubEdges1-offset2-1)];
00465                   Edge *ee=e->getPtr();
00466                   ee->incrRef(); ee->declareOn();
00467                   pushBack(new ElementaryEdge(ee,!(direct1^direction11)));
00468                 }
00469             }
00470         }
00471     }
00472 }
00473 
00474 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
00475 {
00476   int nbOfNodesInPg=0;
00477   bool presenceOfQuadratic=presenceOfQuadraticEdge();
00478   conn.push_back(presenceOfQuadratic?NORM_QPOLYG:NORM_POLYGON);
00479   for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++)
00480     {
00481       Node *tmp=0;
00482       tmp=(*it)->getStartNode();
00483       std::map<INTERP_KERNEL::Node *,int>::const_iterator it1=mapp.find(tmp);
00484       conn.push_back((*it1).second);
00485       nbOfNodesInPg++;
00486     }
00487   if(presenceOfQuadratic)
00488     {
00489       int j=0;
00490       int off=offset+((int)addCoordsQuadratic.size())/2;
00491       for(std::list<ElementaryEdge *>::const_iterator it=_sub_edges.begin();it!=_sub_edges.end();it++,j++,nbOfNodesInPg++)
00492         {
00493           INTERP_KERNEL::Node *node=(*it)->getPtr()->buildRepresentantOfMySelf();
00494           node->unApplySimilarity(xBary,yBary,fact);
00495           addCoordsQuadratic.push_back((*node)[0]);
00496           addCoordsQuadratic.push_back((*node)[1]);
00497           conn.push_back(off+j);
00498           node->decrRef();
00499         }
00500     }
00501   connI.push_back(connI.back()+nbOfNodesInPg+1);
00502 }
00503 
00508 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)
00509 {
00510   double xBaryBB, yBaryBB;
00511   double fact=normalizeExt(&other, xBaryBB, yBaryBB);
00512   //Locate 'this' relative to 'other'
00513   other.performLocatingOperation(*this);
00514   std::vector<QuadraticPolygon *> res=buildIntersectionPolygons(other,*this);
00515   for(std::vector<QuadraticPolygon *>::iterator it=res.begin();it!=res.end();it++)
00516     {
00517       (*it)->appendCrudeData(mapp,xBaryBB,yBaryBB,fact,offset,addCoordsQuadratic,conn,connI);
00518       nbThis.push_back(idThis);
00519       nbOther.push_back(idOther);
00520       delete *it;
00521     }
00522   unApplyGlobalSimilarityExt(other,xBaryBB,yBaryBB,fact);
00523 }
00524 
00529 double QuadraticPolygon::intersectWithAbs1D(QuadraticPolygon& other, bool& isColinear)
00530 {
00531   double ret = 0., xBaryBB, yBaryBB;
00532   double fact = normalize(&other, xBaryBB, yBaryBB);
00533 
00534   QuadraticPolygon cpyOfThis(*this);
00535   QuadraticPolygon cpyOfOther(other);
00536   int nbOfSplits = 0;
00537   SplitPolygonsEachOther(cpyOfThis, cpyOfOther, nbOfSplits);
00538   //At this point cpyOfThis and cpyOfOther have been splited at maximum edge so that in/out can been done.
00539   performLocatingOperation(cpyOfOther);
00540   isColinear = false;
00541   for(std::list<ElementaryEdge *>::const_iterator it=cpyOfOther._sub_edges.begin();it!=cpyOfOther._sub_edges.end();it++)
00542     {
00543       switch((*it)->getLoc())
00544         {
00545         case FULL_IN_1:
00546           {
00547             ret += fabs((*it)->getPtr()->getCurveLength());
00548             break;
00549           }
00550         case FULL_ON_1:
00551           {
00552             isColinear=true;
00553             ret += fabs((*it)->getPtr()->getCurveLength());
00554             break;
00555           }
00556         default:
00557           {
00558           }
00559         }
00560     }
00561   return ret * fact;
00562 }
00563 
00567 double QuadraticPolygon::intersectWithAbs(QuadraticPolygon& other, double* barycenter)
00568 {
00569   double ret=0.,bary[2],area,xBaryBB,yBaryBB;
00570   barycenter[0] = barycenter[1] = 0.;
00571   double fact=normalize(&other,xBaryBB,yBaryBB);
00572   std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
00573   for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
00574     {
00575       area=fabs((*iter)->getArea());
00576       (*iter)->getBarycenter(bary);
00577       delete *iter;
00578       ret+=area;
00579       barycenter[0] += bary[0]*area;
00580       barycenter[1] += bary[1]*area;
00581     }
00582   if ( ret > std::numeric_limits<double>::min() )
00583     {
00584       barycenter[0]=barycenter[0]/ret*fact+xBaryBB;
00585       barycenter[1]=barycenter[1]/ret*fact+yBaryBB;
00586       
00587     }
00588   return ret*fact*fact;
00589 }
00590 
00596 double QuadraticPolygon::intersectWith(const QuadraticPolygon& other) const
00597 {
00598   double ret=0.;
00599   std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
00600   for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
00601     {
00602       ret+=fabs((*iter)->getArea());
00603       delete *iter;
00604     }
00605   return ret;
00606 }
00607 
00613 double QuadraticPolygon::intersectWith(const QuadraticPolygon& other, double* barycenter) const
00614 {
00615   double ret=0., bary[2];
00616   barycenter[0] = barycenter[1] = 0.;
00617   std::vector<QuadraticPolygon *> polygs=intersectMySelfWith(other);
00618   for(std::vector<QuadraticPolygon *>::iterator iter=polygs.begin();iter!=polygs.end();iter++)
00619     {
00620       double area = fabs((*iter)->getArea());
00621       (*iter)->getBarycenter(bary);
00622       delete *iter;
00623       ret+=area;
00624       barycenter[0] += bary[0]*area;
00625       barycenter[1] += bary[1]*area;
00626     }
00627   if ( ret > std::numeric_limits<double>::min() )
00628     {
00629       barycenter[0] /= ret;
00630       barycenter[1] /= ret;
00631     }
00632   return ret;
00633 }
00634 
00640 void QuadraticPolygon::intersectForPerimeter(const QuadraticPolygon& other, double& perimeterThisPart, double& perimeterOtherPart, double& perimeterCommonPart) const
00641 {
00642   perimeterThisPart=0.; perimeterOtherPart=0.; perimeterCommonPart=0.;
00643   QuadraticPolygon cpyOfThis(*this);
00644   QuadraticPolygon cpyOfOther(other); int nbOfSplits=0;
00645   SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits);
00646   performLocatingOperation(cpyOfOther);
00647   other.performLocatingOperation(cpyOfThis);
00648   cpyOfThis.dispatchPerimeterExcl(perimeterThisPart,perimeterCommonPart);
00649   cpyOfOther.dispatchPerimeterExcl(perimeterOtherPart,perimeterCommonPart);
00650   perimeterCommonPart/=2.;
00651 }
00652 
00663 void QuadraticPolygon::intersectForPerimeterAdvanced(const QuadraticPolygon& other, std::vector< double >& polThis, std::vector< double >& polOther) const
00664 {
00665   polThis.resize(size());
00666   polOther.resize(other.size());
00667   IteratorOnComposedEdge it1(const_cast<QuadraticPolygon *>(this));
00668   int edgeId=0;
00669   for(it1.first();!it1.finished();it1.next(),edgeId++)
00670     {
00671       ElementaryEdge* curE1=it1.current();
00672       QuadraticPolygon cpyOfOther(other);
00673       QuadraticPolygon tmp;
00674       tmp.pushBack(curE1->clone());
00675       int tmp2;
00676       SplitPolygonsEachOther(tmp,cpyOfOther,tmp2);
00677       other.performLocatingOperation(tmp);
00678       tmp.dispatchPerimeter(polThis[edgeId]);
00679     }
00680   //
00681   IteratorOnComposedEdge it2(const_cast<QuadraticPolygon *>(&other));
00682   edgeId=0;
00683   for(it2.first();!it2.finished();it2.next(),edgeId++)
00684     {
00685       ElementaryEdge* curE2=it2.current();
00686       QuadraticPolygon cpyOfThis(*this);
00687       QuadraticPolygon tmp;
00688       tmp.pushBack(curE2->clone());
00689       int tmp2;
00690       SplitPolygonsEachOther(tmp,cpyOfThis,tmp2);
00691       performLocatingOperation(tmp);
00692       tmp.dispatchPerimeter(polOther[edgeId]);
00693     }
00694 }
00695 
00696 
00703 void QuadraticPolygon::intersectForPoint(const QuadraticPolygon& other, std::vector< int >& numberOfCreatedPointsPerEdge) const
00704 {
00705   numberOfCreatedPointsPerEdge.resize(size());
00706   IteratorOnComposedEdge it1(const_cast<QuadraticPolygon *>(this));
00707   int edgeId=0;
00708   for(it1.first();!it1.finished();it1.next(),edgeId++)
00709     {
00710       ElementaryEdge* curE1=it1.current();
00711       QuadraticPolygon cpyOfOther(other);
00712       QuadraticPolygon tmp;
00713       tmp.pushBack(curE1->clone());
00714       int tmp2;
00715       SplitPolygonsEachOther(tmp,cpyOfOther,tmp2);
00716       numberOfCreatedPointsPerEdge[edgeId]=tmp.recursiveSize()-1;
00717     }
00718 }
00719 
00725 std::vector<QuadraticPolygon *> QuadraticPolygon::intersectMySelfWith(const QuadraticPolygon& other) const
00726 {
00727   QuadraticPolygon cpyOfThis(*this);
00728   QuadraticPolygon cpyOfOther(other); int nbOfSplits=0;
00729   SplitPolygonsEachOther(cpyOfThis,cpyOfOther,nbOfSplits);
00730   //At this point cpyOfThis and cpyOfOther have been splited at maximum edge so that in/out can been done.
00731   performLocatingOperation(cpyOfOther);
00732   return other.buildIntersectionPolygons(cpyOfThis,cpyOfOther);
00733 }
00734 
00740 void QuadraticPolygon::SplitPolygonsEachOther(QuadraticPolygon& pol1, QuadraticPolygon& pol2, int& nbOfSplits)
00741 {
00742   IteratorOnComposedEdge it1(&pol1),it2(&pol2);
00743   MergePoints merge;
00744   ComposedEdge *c1=new ComposedEdge;
00745   ComposedEdge *c2=new ComposedEdge;
00746   for(it2.first();!it2.finished();it2.next())
00747     {
00748       ElementaryEdge* curE2=it2.current();
00749       if(!curE2->isThereStartPoint())
00750         it1.first();
00751       else
00752         it1=curE2->getIterator();
00753       for(;!it1.finished();)
00754         {
00755           
00756           ElementaryEdge* curE1=it1.current();
00757           merge.clear(); nbOfSplits++;
00758           if(curE1->getPtr()->intersectWith(curE2->getPtr(),merge,*c1,*c2))
00759             {
00760               if(!curE1->getDirection()) c1->reverse();
00761               if(!curE2->getDirection()) c2->reverse();
00762               UpdateNeighbours(merge,it1,it2,c1,c2);
00763               //Substitution of simple edge by sub-edges.
00764               delete curE1; // <-- destroying simple edge coming from pol1
00765               delete curE2; // <-- destroying simple edge coming from pol2
00766               it1.insertElemEdges(c1,true);// <-- 2nd param is true to go next.
00767               it2.insertElemEdges(c2,false);// <-- 2nd param is false to avoid to go next.
00768               curE2=it2.current();
00769               //
00770               it1.assignMySelfToAllElems(c2);//To avoid that others
00771               SoftDelete(c1);
00772               SoftDelete(c2);
00773               c1=new ComposedEdge;
00774               c2=new ComposedEdge;
00775             }
00776           else
00777             {
00778               UpdateNeighbours(merge,it1,it2,curE1,curE2);
00779               it1.next();
00780             }
00781         }
00782     }
00783   Delete(c1);
00784   Delete(c2);
00785 }
00786 
00787 void QuadraticPolygon::performLocatingOperation(QuadraticPolygon& pol2) const
00788 {
00789   IteratorOnComposedEdge it(&pol2);
00790   TypeOfEdgeLocInPolygon loc=FULL_ON_1;
00791   for(it.first();!it.finished();it.next())
00792     {
00793       ElementaryEdge *cur=it.current();
00794       loc=cur->locateFullyMySelf(*this,loc);
00795     }
00796 }
00797 
00805 std::vector<QuadraticPolygon *> QuadraticPolygon::buildIntersectionPolygons(const QuadraticPolygon& pol1, const QuadraticPolygon& pol2) const
00806 {
00807   std::vector<QuadraticPolygon *> ret;
00808   std::list<QuadraticPolygon *> pol2Zip=pol2.zipConsecutiveInSegments();
00809   if(!pol2Zip.empty())
00810     closePolygons(pol2Zip,pol1,ret);
00811   else
00812     {//borders of pol2 do not cross pol1,and pol2 borders are outside of pol1. That is to say, either pol2 and pol1
00813       //do not overlap or  pol1 is fully inside pol2. So in the first case no intersection, in the other case
00814       //the intersection is pol1.
00815       ElementaryEdge *e1FromPol1=pol1[0];
00816       TypeOfEdgeLocInPolygon loc=FULL_ON_1;
00817       loc=e1FromPol1->locateFullyMySelf(*this,loc);
00818       if(loc==FULL_IN_1)
00819         ret.push_back(new QuadraticPolygon(pol1));
00820     }
00821   return ret;
00822 }
00823 
00828 std::list<QuadraticPolygon *> QuadraticPolygon::zipConsecutiveInSegments() const
00829 {
00830   std::list<QuadraticPolygon *> ret;
00831   IteratorOnComposedEdge it((ComposedEdge *)this);
00832   int nbOfTurns=recursiveSize();
00833   int i=0;
00834   if(!it.goToNextInOn(false,i,nbOfTurns))
00835     return ret;
00836   i=0;
00837   //
00838   while(i<nbOfTurns)
00839     {
00840       QuadraticPolygon *tmp1=new QuadraticPolygon;
00841       TypeOfEdgeLocInPolygon loc=it.current()->getLoc();
00842       while(loc!=FULL_OUT_1 && i<nbOfTurns)
00843         {
00844           ElementaryEdge *tmp3=it.current()->clone();
00845           tmp1->pushBack(tmp3);
00846           it.nextLoop(); i++;
00847           loc=it.current()->getLoc();
00848         }
00849       if(tmp1->empty())
00850         {
00851           delete tmp1;
00852           continue;
00853         }
00854       ret.push_back(tmp1);
00855       it.goToNextInOn(true,i,nbOfTurns);
00856     }
00857   return ret;
00858 }
00859 
00866 void QuadraticPolygon::closePolygons(std::list<QuadraticPolygon *>& pol2Zip, const QuadraticPolygon& pol1,
00867                                      std::vector<QuadraticPolygon *>& results) const
00868 {
00869   bool directionKnownInPol1=false;
00870   bool directionInPol1;
00871   for(std::list<QuadraticPolygon *>::iterator iter=pol2Zip.begin();iter!=pol2Zip.end();)
00872     {
00873       if((*iter)->completed())
00874         {
00875           results.push_back(*iter);
00876           directionKnownInPol1=false;
00877           iter=pol2Zip.erase(iter);
00878           continue;
00879         }
00880       if(!directionKnownInPol1)
00881         {
00882           if(!(*iter)->amIAChanceToBeCompletedBy(pol1,*this,directionInPol1))
00883             { delete *iter; iter=pol2Zip.erase(iter); continue; }
00884           else
00885             directionKnownInPol1=true;
00886         }
00887       std::list<QuadraticPolygon *>::iterator iter2=iter; iter2++;
00888       std::list<QuadraticPolygon *>::iterator iter3=(*iter)->fillAsMuchAsPossibleWith(pol1,iter2,pol2Zip.end(),directionInPol1);
00889       if(iter3!=pol2Zip.end())
00890         {
00891           (*iter)->pushBack(*iter3);
00892           SoftDelete(*iter3);
00893           pol2Zip.erase(iter3);
00894         }
00895     }
00896 }
00897 
00901 bool QuadraticPolygon::amIAChanceToBeCompletedBy(const QuadraticPolygon& pol1Splitted,const QuadraticPolygon& pol2NotSplitted, bool& direction)
00902 {
00903   IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(&pol1Splitted));
00904   bool found=false;
00905   Node *n=getEndNode();
00906   ElementaryEdge *cur=it.current();
00907   for(it.first();!it.finished() && !found;)
00908     {
00909       cur=it.current();
00910       found=(cur->getStartNode()==n);
00911       if(!found)
00912         it.next();
00913     }
00914   if(!found)
00915     throw Exception("Internal error : polygons uncompatible each others. Should never happend");
00916   //Ok we found correspondance between this and pol1. Searching for right direction to close polygon.
00917   ElementaryEdge *e=_sub_edges.back();
00918   if(e->getLoc()==FULL_ON_1)
00919     {
00920       if(e->getPtr()==cur->getPtr())
00921         {
00922           direction=false;
00923           it.previousLoop();
00924           cur=it.current();
00925           Node *repr=cur->getPtr()->buildRepresentantOfMySelf();
00926           bool ret=pol2NotSplitted.isInOrOut(repr);
00927           repr->decrRef();
00928           return ret;
00929         }
00930       else
00931         {
00932           direction=true;
00933           Node *repr=cur->getPtr()->buildRepresentantOfMySelf();
00934           bool ret=pol2NotSplitted.isInOrOut(repr);
00935           repr->decrRef();
00936           return ret;
00937         }
00938     }
00939   else
00940     direction=cur->locateFullyMySelfAbsolute(pol2NotSplitted)==FULL_IN_1;
00941   return true;
00942 }
00943 
00947 std::list<QuadraticPolygon *>::iterator QuadraticPolygon::fillAsMuchAsPossibleWith(const QuadraticPolygon& pol1Splitted,
00948                                                                                    std::list<QuadraticPolygon *>::iterator iStart,
00949                                                                                    std::list<QuadraticPolygon *>::iterator iEnd,
00950                                                                                    bool direction)
00951 {
00952   IteratorOnComposedEdge it(const_cast<QuadraticPolygon *>(&pol1Splitted));
00953   bool found=false;
00954   Node *n=getEndNode();
00955   ElementaryEdge *cur;
00956   for(it.first();!it.finished() && !found;)
00957     {
00958       cur=it.current();
00959       found=(cur->getStartNode()==n);
00960       if(!found)
00961         it.next();
00962     }
00963   if(!direction)
00964     it.previousLoop();
00965   Node *nodeToTest;
00966   std::list<QuadraticPolygon *>::iterator ret;
00967   do
00968     {
00969       cur=it.current();
00970       ElementaryEdge *tmp=cur->clone();
00971       if(!direction)
00972         tmp->reverse();
00973       pushBack(tmp);
00974       nodeToTest=tmp->getEndNode();
00975       direction?it.nextLoop():it.previousLoop();
00976       ret=CheckInList(nodeToTest,iStart,iEnd);
00977       if(completed())
00978         return iEnd;
00979     }
00980   while(ret==iEnd);
00981   return ret;
00982 }
00983 
00984 std::list<QuadraticPolygon *>::iterator QuadraticPolygon::CheckInList(Node *n, std::list<QuadraticPolygon *>::iterator iStart,
00985                                                                       std::list<QuadraticPolygon *>::iterator iEnd)
00986 {
00987   for(std::list<QuadraticPolygon *>::iterator iter=iStart;iter!=iEnd;iter++)
00988     if((*iter)->isNodeIn(n))
00989       return iter;
00990   return iEnd;
00991 }