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