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