0032328: Missing include of TopoDS_Edge.hxx in ShapeUpgrade_UnifySameDomain.hxx
[occt.git] / src / ShapeUpgrade / ShapeUpgrade_UnifySameDomain.cxx
1 // Created on: 2012-06-09
2 // Created by: jgv@ROLEX
3 // Copyright (c) 2012-2014 OPEN CASCADE SAS
4 //
5 // This file is part of Open CASCADE Technology software library.
6 //
7 // This library is free software; you can redistribute it and/or modify it under
8 // the terms of the GNU Lesser General Public License version 2.1 as published
9 // by the Free Software Foundation, with special exception defined in the file
10 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
11 // distribution for complete text of the license and disclaimer of any warranty.
12 //
13 // Alternatively, this file may be used under the terms of Open CASCADE
14 // commercial license or contractual agreement.
15
16 #include <ShapeUpgrade_UnifySameDomain.hxx>
17
18 #include <BRep_Builder.hxx>
19 #include <BRep_Tool.hxx>
20 #include <BRepLib.hxx>
21 #include <BRepLib_MakeEdge.hxx>
22 #include <BRepLib_MakeFace.hxx>
23 #include <BRepTopAdaptor_TopolTool.hxx>
24 #include <GC_MakeCircle.hxx>
25 #include <Geom2d_Line.hxx>
26 #include <GCE2d_MakeLine.hxx>
27 #include <Geom2d_TrimmedCurve.hxx>
28 #include <Geom2dConvert.hxx>
29 #include <Geom2dConvert_CompCurveToBSplineCurve.hxx>
30 #include <Geom_BezierCurve.hxx>
31 #include <Geom_BSplineCurve.hxx>
32 #include <Geom_Circle.hxx>
33 #include <Geom_Ellipse.hxx>
34 #include <Geom_Hyperbola.hxx>
35 #include <Geom_Parabola.hxx>
36 #include <Geom_CylindricalSurface.hxx>
37 #include <Geom_ElementarySurface.hxx>
38 #include <Geom_Line.hxx>
39 #include <Geom_OffsetSurface.hxx>
40 #include <Geom_Plane.hxx>
41 #include <Geom_RectangularTrimmedSurface.hxx>
42 #include <Geom_Surface.hxx>
43 #include <Geom_SurfaceOfLinearExtrusion.hxx>
44 #include <Geom_SurfaceOfRevolution.hxx>
45 #include <Geom_SweptSurface.hxx>
46 #include <Geom_TrimmedCurve.hxx>
47 #include <GeomAdaptor_Surface.hxx>
48 #include <GeomConvert.hxx>
49 #include <GeomConvert_CompCurveToBSplineCurve.hxx>
50 #include <GeomLib_IsPlanarSurface.hxx>
51 #include <gp_Cylinder.hxx>
52 #include <gp_Dir.hxx>
53 #include <gp_Lin.hxx>
54 #include <IntPatch_ImpImpIntersection.hxx>
55 #include <ShapeAnalysis_Edge.hxx>
56 #include <ShapeAnalysis_WireOrder.hxx>
57 #include <ShapeAnalysis_Surface.hxx>
58 #include <ShapeBuild_Edge.hxx>
59 #include <ShapeBuild_ReShape.hxx>
60 #include <ShapeConstruct_ProjectCurveOnSurface.hxx>
61 #include <ShapeFix_Edge.hxx>
62 #include <ShapeFix_Face.hxx>
63 #include <ShapeFix_Shell.hxx>
64 #include <ShapeFix_Wire.hxx>
65 #include <Standard_Type.hxx>
66 #include <TColGeom2d_Array1OfBSplineCurve.hxx>
67 #include <TColGeom2d_HArray1OfBSplineCurve.hxx>
68 #include <TColGeom2d_SequenceOfBoundedCurve.hxx>
69 #include <TColGeom2d_SequenceOfCurve.hxx>
70 #include <TColGeom_Array1OfBSplineCurve.hxx>
71 #include <TColGeom_HArray1OfBSplineCurve.hxx>
72 #include <TColGeom_SequenceOfSurface.hxx>
73 #include <TColStd_Array1OfReal.hxx>
74 #include <TColStd_MapOfInteger.hxx>
75 #include <TColStd_SequenceOfBoolean.hxx>
76 #include <TopExp.hxx>
77 #include <TopExp_Explorer.hxx>
78 #include <TopoDS.hxx>
79 #include <TopoDS_Edge.hxx>
80 #include <TopoDS_Face.hxx>
81 #include <TopoDS_Shape.hxx>
82 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
83 #include <TopTools_IndexedMapOfShape.hxx>
84 #include <TopTools_ListIteratorOfListOfShape.hxx>
85 #include <TopTools_MapOfShape.hxx>
86 #include <TopTools_SequenceOfShape.hxx>
87 #include <gp_Circ.hxx>
88 #include <BRepAdaptor_Curve.hxx>
89 #include <BRepAdaptor_Curve2d.hxx>
90 #include <BRepAdaptor_Surface.hxx>
91 #include <gp_Vec2d.hxx>
92 #include <Extrema_ExtPS.hxx>
93 #include <BRepTools.hxx>
94 #include <BRepTopAdaptor_FClass2d.hxx>
95 #include <ElCLib.hxx>
96 #include <BRepBuilderAPI_MakeEdge.hxx>
97 #include <BRepBuilderAPI_MakeWire.hxx>
98 #include <BRepBuilderAPI_MakeFace.hxx>
99 #include <GCPnts_AbscissaPoint.hxx>
100 #include <ElSLib.hxx>
101 #include <GeomProjLib.hxx>
102
103 IMPLEMENT_STANDARD_RTTIEXT(ShapeUpgrade_UnifySameDomain,Standard_Transient)
104
105 static void SplitWire (const TopoDS_Wire&                theWire,
106                        const TopoDS_Face&                theFace,
107                        const TopTools_IndexedMapOfShape& theVmap,
108                        TopTools_SequenceOfShape&         theWireSeq);
109
110 static Standard_Real TrueValueOfOffset(const Standard_Real theValue,
111                                        const Standard_Real thePeriod)
112 {
113   if (theValue > 0.)
114     return thePeriod;
115
116   return (-thePeriod);
117 }
118
119 //=======================================================================
120 //function : UpdateBoundaries
121 //purpose  : auxiliary
122 //=======================================================================
123 static void UpdateBoundaries(const Handle(Geom2d_Curve)& thePCurve,
124                              const Standard_Real         theFirst,
125                              const Standard_Real         theLast,
126                              Standard_Real& theUfirst,
127                              Standard_Real& theUlast)
128 {
129   const Standard_Integer NbSamples = 4;
130   Standard_Real Delta = (theLast - theFirst)/NbSamples;
131
132   for (Standard_Integer i = 0; i <= NbSamples; i++)
133   {
134     Standard_Real aParam = theFirst + i*Delta;
135     gp_Pnt2d aPoint = thePCurve->Value(aParam);
136     if (aPoint.X() < theUfirst)
137       theUfirst = aPoint.X();
138     if (aPoint.X() > theUlast)
139       theUlast  = aPoint.X();
140   }
141 }
142
143 static Standard_Boolean TryMakeLine(const Handle(Geom2d_Curve)& thePCurve,
144                                     const Standard_Real         theFirst,
145                                     const Standard_Real         theLast,
146                                     Handle(Geom2d_Line)&        theLine)
147 {
148   gp_Pnt2d aFirstPnt = thePCurve->Value (theFirst);
149   gp_Pnt2d aLastPnt  = thePCurve->Value (theLast);
150   gp_Vec2d aVec (aFirstPnt, aLastPnt);
151   Standard_Real aSqLen = aVec.SquareMagnitude();
152   Standard_Real aSqParamLen = (theLast - theFirst)*(theLast - theFirst);
153   if (Abs(aSqLen - aSqParamLen) > Precision::Confusion())
154     return Standard_False;
155
156   gp_Dir2d aDir = aVec;
157   gp_Vec2d anOffset = -aDir;
158   anOffset *= theFirst;
159   gp_Pnt2d anOrigin = aFirstPnt.Translated(anOffset);
160   gp_Lin2d aLin (anOrigin, aDir);
161
162   const Standard_Integer NbSamples = 10;
163   Standard_Real aDelta = (theLast - theFirst)/NbSamples;
164   for (Standard_Integer i = 1; i < NbSamples; i++)
165   {
166     Standard_Real aParam = theFirst + i*aDelta;
167     gp_Pnt2d aPnt = thePCurve->Value(aParam);
168     Standard_Real aDist = aLin.Distance (aPnt);
169     if (aDist > Precision::Confusion())
170       return Standard_False;
171   }
172
173   theLine = new Geom2d_Line (aLin);
174   return Standard_True;
175 }
176
177 static void RemoveEdgeFromMap(const TopoDS_Edge& theEdge,
178                               TopTools_IndexedDataMapOfShapeListOfShape& theVEmap)
179 {
180   TopoDS_Vertex VV [2];
181   TopExp::Vertices(theEdge, VV[0], VV[1]);
182   for (Standard_Integer i = 0; i < 2; i++)
183   {
184     TopTools_ListOfShape& Elist = theVEmap.ChangeFromKey(VV[i]);
185     TopTools_ListIteratorOfListOfShape itl(Elist);
186     while (itl.More())
187     {
188       const TopoDS_Shape& anEdge = itl.Value();
189       if (anEdge.IsSame(theEdge))
190         Elist.Remove(itl);
191       else
192         itl.Next();
193     }
194   }
195 }
196
197 static Standard_Real ComputeMinEdgeSize(const TopTools_SequenceOfShape& theEdges,
198                                         const TopoDS_Face&              theRefFace,
199                                         TopTools_MapOfShape&            theEdgesMap)
200 {
201   Standard_Real MinSize = RealLast();
202   
203   for (Standard_Integer ind = 1; ind <= theEdges.Length(); ind++)
204   {
205     const TopoDS_Edge& anEdge = TopoDS::Edge(theEdges(ind));
206     theEdgesMap.Add(anEdge);
207     TopoDS_Vertex V1, V2;
208     TopExp::Vertices(anEdge, V1, V2);
209     BRepAdaptor_Curve2d BAcurve2d(anEdge, theRefFace);
210     gp_Pnt2d FirstP2d = BAcurve2d.Value(BAcurve2d.FirstParameter());
211     gp_Pnt2d LastP2d  = BAcurve2d.Value(BAcurve2d.LastParameter());
212     Standard_Real aSqDist;
213     if (V1.IsSame(V2) &&
214         !BRep_Tool::Degenerated(anEdge))
215     {
216       gp_Pnt2d MidP2d = BAcurve2d.Value((BAcurve2d.FirstParameter()+BAcurve2d.LastParameter())/2);
217       aSqDist = FirstP2d.SquareDistance(MidP2d);
218     }
219     else
220       aSqDist = FirstP2d.SquareDistance(LastP2d);
221     if (aSqDist < MinSize)
222       MinSize = aSqDist;
223   }
224   MinSize = Sqrt(MinSize);
225   return MinSize;
226 }
227
228 static void FindFaceWithMaxUbound(const TopTools_SequenceOfShape& theFaces,
229                                   const TopoDS_Face&              theRefFace,
230                                   const TopTools_MapOfShape&      theEdgesMap,
231                                   Standard_Real&                  theUmax,
232                                   Standard_Integer&               theIndFaceMax)
233 {
234   theUmax = RealFirst();
235   theIndFaceMax = 0;
236   for (Standard_Integer ii = 1; ii <= theFaces.Length(); ii++)
237   {
238     const TopoDS_Face& face_ii = TopoDS::Face(theFaces(ii));
239     Standard_Real Ufirst = RealLast(), Ulast = RealFirst();
240     TopExp_Explorer Explo(face_ii, TopAbs_EDGE);
241     for (; Explo.More(); Explo.Next())
242     {
243       const TopoDS_Edge& anEdge = TopoDS::Edge(Explo.Current());
244       if (!theEdgesMap.Contains(anEdge))
245         continue;
246       Standard_Real fpar, lpar;
247       Handle(Geom2d_Curve) aPCurve = BRep_Tool::CurveOnSurface(anEdge, theRefFace, fpar, lpar);
248       UpdateBoundaries(aPCurve, fpar, lpar, Ufirst, Ulast);
249     }
250     if (Ulast > theUmax)
251     {
252       theUmax = Ulast;
253       theIndFaceMax = ii;
254     }
255   }
256 }
257
258 static void RelocatePCurvesToNewUorigin(const TopTools_SequenceOfShape& theEdges,
259                                         const TopoDS_Shape&             theFirstFace,
260                                         const TopoDS_Face&              theRefFace,
261                                         const Standard_Real             theCoordTol,
262                                         const Standard_Real             theUperiod,
263                                         TopTools_IndexedDataMapOfShapeListOfShape& theVEmap,
264                                         NCollection_DataMap<TopoDS_Shape, Handle(Geom2d_Curve)>& theEdgeNewPCurve,
265                                         TopTools_MapOfShape&      theUsedEdges)
266 {
267   TopTools_MapOfShape EdgesOfFirstFace;
268   TopExp::MapShapes(theFirstFace, EdgesOfFirstFace);
269   
270   for (;;) //walk by contours
271   {
272     //try to find the start edge that:
273     //1. belongs to outer edges of first face;
274     //2. has minimum U-coord of its start point
275     TopoDS_Edge StartEdge;
276     TopAbs_Orientation anOr = TopAbs_FORWARD;
277     Standard_Real anUmin = RealLast();
278     for (Standard_Integer ii = 1; ii <= theEdges.Length(); ii++)
279     {
280       const TopoDS_Edge& anEdge = TopoDS::Edge(theEdges(ii));
281       if (theUsedEdges.Contains(anEdge))
282         continue;
283       if (EdgesOfFirstFace.Contains(anEdge))
284       {
285         if (StartEdge.IsNull())
286         {
287           StartEdge = anEdge;
288           BRepAdaptor_Curve2d StartBAcurve(StartEdge, theRefFace);
289           Standard_Real aFirstParam, aLastParam;
290           if (StartEdge.Orientation() == TopAbs_FORWARD)
291           {
292             aFirstParam = StartBAcurve.FirstParameter();
293             aLastParam  = StartBAcurve.LastParameter();
294           }
295           else
296           {
297             aFirstParam = StartBAcurve.LastParameter();
298             aLastParam  = StartBAcurve.FirstParameter();
299           }
300           gp_Pnt2d aFirstPoint = StartBAcurve.Value(aFirstParam);
301           gp_Pnt2d aLastPoint  = StartBAcurve.Value(aLastParam);
302           if (aFirstPoint.X() < aLastPoint.X())
303           {
304             anUmin = aFirstPoint.X();
305             anOr = TopAbs_FORWARD;
306           }
307           else
308           {
309             anUmin = aLastPoint.X();
310             anOr = TopAbs_REVERSED;
311           }
312         }
313         else
314         {
315           BRepAdaptor_Curve2d aBAcurve(anEdge, theRefFace);
316           Standard_Real aFirstParam, aLastParam;
317           if (anEdge.Orientation() == TopAbs_FORWARD)
318           {
319             aFirstParam = aBAcurve.FirstParameter();
320             aLastParam  = aBAcurve.LastParameter();
321           }
322           else
323           {
324             aFirstParam = aBAcurve.LastParameter();
325             aLastParam  = aBAcurve.FirstParameter();
326           }
327           gp_Pnt2d aFirstPoint = aBAcurve.Value(aFirstParam);
328           gp_Pnt2d aLastPoint  = aBAcurve.Value(aLastParam);
329           if (aFirstPoint.X() < anUmin)
330           {
331             StartEdge = anEdge;
332             anUmin = aFirstPoint.X();
333             anOr = TopAbs_FORWARD;
334           }
335           if (aLastPoint.X() < anUmin)
336           {
337             StartEdge = anEdge;
338             anUmin = aLastPoint.X();
339             anOr = TopAbs_REVERSED;
340           }
341         }
342       } //if (EdgesOfFirstFace.Contains(anEdge))
343     } //for (Standard_Integer ii = 1; ii <= edges.Length(); ii++)
344     
345     if (StartEdge.IsNull()) //all contours are passed
346       break;
347     
348     TopoDS_Vertex StartVertex = (anOr == TopAbs_FORWARD)?
349       TopExp::FirstVertex(StartEdge, Standard_True) : TopExp::LastVertex(StartEdge, Standard_True);
350     TopoDS_Edge CurEdge = StartEdge;
351     Standard_Real fpar, lpar;
352     Handle(Geom2d_Curve) CurPCurve = BRep_Tool::CurveOnSurface(CurEdge, theRefFace, fpar, lpar);
353     CurPCurve = new Geom2d_TrimmedCurve(CurPCurve, fpar, lpar);
354     theEdgeNewPCurve.Bind(CurEdge, CurPCurve);
355     if (CurEdge.Orientation() == TopAbs_REVERSED)
356     { Standard_Real tmp = fpar; fpar = lpar; lpar = tmp; }
357     //Standard_Real StartParam = (anOr == TopAbs_FORWARD)? fpar : lpar;
358     Standard_Real CurParam = (anOr == TopAbs_FORWARD)? lpar : fpar;
359     //gp_Pnt2d StartPoint = CurPCurve->Value(StartParam);
360     gp_Pnt2d CurPoint = CurPCurve->Value(CurParam);
361     for (;;) //collect pcurves of a contour
362     {
363       RemoveEdgeFromMap(CurEdge, theVEmap);
364       theUsedEdges.Add(CurEdge);
365       TopoDS_Vertex CurVertex = (anOr == TopAbs_FORWARD)?
366         TopExp::LastVertex(CurEdge, Standard_True) : TopExp::FirstVertex(CurEdge, Standard_True);
367       
368       const TopTools_ListOfShape& Elist = theVEmap.FindFromKey(CurVertex);
369       if (Elist.IsEmpty())
370         break; //end of contour in 3d
371       
372       TopTools_ListIteratorOfListOfShape itl(Elist);
373       for (; itl.More(); itl.Next())
374       {
375         const TopoDS_Edge& anEdge = TopoDS::Edge(itl.Value());
376         if (anEdge.IsSame(CurEdge))
377           continue;
378         TopoDS_Vertex aFirstVertex = (anOr == TopAbs_FORWARD)?
379           TopExp::FirstVertex(anEdge, Standard_True) : TopExp::LastVertex(anEdge, Standard_True);
380         if (!aFirstVertex.IsSame(CurVertex)) //may be if CurVertex is deg.vertex
381           continue;
382         
383         Handle(Geom2d_Curve) aPCurve = BRep_Tool::CurveOnSurface(anEdge, theRefFace, fpar, lpar);
384         aPCurve = new Geom2d_TrimmedCurve(aPCurve, fpar, lpar);
385         if (anEdge.Orientation() == TopAbs_REVERSED)
386         { Standard_Real tmp = fpar; fpar = lpar; lpar = tmp; }
387         Standard_Real aParam = (anOr == TopAbs_FORWARD)? fpar : lpar;
388         gp_Pnt2d aPoint = aPCurve->Value(aParam);
389         Standard_Real anOffset = CurPoint.X() - aPoint.X();
390         if (!(Abs(anOffset) < theCoordTol ||
391               Abs(Abs(anOffset) - theUperiod) < theCoordTol))
392           continue; //may be if CurVertex is deg.vertex
393         
394         if (Abs(anOffset) > theUperiod/2)
395         {
396           anOffset = TrueValueOfOffset(anOffset, theUperiod);
397           gp_Vec2d aVec(anOffset, 0.);
398           Handle(Geom2d_Curve) aNewPCurve = Handle(Geom2d_Curve)::DownCast(aPCurve->Copy());
399           aNewPCurve->Translate(aVec);
400           aPCurve = aNewPCurve;
401         }
402         theEdgeNewPCurve.Bind(anEdge, aPCurve);
403         CurEdge = anEdge;
404         TopAbs_Orientation CurOr = TopAbs::Compose(anOr, CurEdge.Orientation());
405         CurParam = (CurOr == TopAbs_FORWARD)?
406           aPCurve->LastParameter() : aPCurve->FirstParameter();
407         CurPoint = aPCurve->Value(CurParam);
408         break;
409       }
410     } //for (;;) (collect pcurves of a contour)
411   } //for (;;) (walk by contours)
412 }
413
414 static void InsertWiresIntoFaces(const TopTools_SequenceOfShape& theWires,
415                                  const TopTools_SequenceOfShape& theFaces,
416                                  const TopoDS_Face&              theRefFace)
417 {
418   BRep_Builder BB;
419   for (Standard_Integer ii = 1; ii <= theWires.Length(); ii++)
420   {
421     const TopoDS_Wire& aWire = TopoDS::Wire(theWires(ii));
422     TopoDS_Iterator iter(aWire);
423     const TopoDS_Edge& anEdge = TopoDS::Edge(iter.Value());
424     BRepAdaptor_Curve2d BAcurve2d(anEdge, theRefFace);
425     gp_Pnt2d aPnt2d = BAcurve2d.Value((BAcurve2d.FirstParameter() + BAcurve2d.LastParameter())/2.);
426     TopoDS_Shape RequiredFace;
427     for (Standard_Integer jj = 1; jj <= theFaces.Length(); jj++)
428     {
429       const TopoDS_Face& aFace = TopoDS::Face(theFaces(jj));
430       BRepTopAdaptor_FClass2d Classifier(aFace, Precision::Confusion());
431       TopAbs_State aStatus = Classifier.Perform(aPnt2d);
432       if (aStatus == TopAbs_IN)
433       {
434         RequiredFace = aFace.Oriented (TopAbs_FORWARD);
435         break;
436       }
437     }
438     if (!RequiredFace.IsNull())
439     {
440       BB.Add(RequiredFace, aWire);
441     }
442     else
443     {
444       Standard_ASSERT_INVOKE ("ShapeUpgrade_UnifySameDomain: wire remains unclassified");
445     }
446   }
447 }
448
449 static TopoDS_Face FindCommonFace(const TopoDS_Edge&  theEdge1,
450                                   const TopoDS_Edge&  theEdge2,
451                                   const TopTools_IndexedDataMapOfShapeListOfShape& theVFmap,
452                                   TopAbs_Orientation& theOrOfE1OnFace,
453                                   TopAbs_Orientation& theOrOfE2OnFace)
454 {
455   TopoDS_Vertex aVertex;
456   TopExp::CommonVertex(theEdge1, theEdge2, aVertex);
457   const TopTools_ListOfShape& Flist = theVFmap.FindFromKey(aVertex);
458   TopTools_ListIteratorOfListOfShape itl(Flist);
459   for (; itl.More(); itl.Next())
460   {
461     TopoDS_Face aFace = TopoDS::Face(itl.Value());
462     aFace.Orientation(TopAbs_FORWARD);
463     Standard_Boolean e1found = Standard_False, e2found = Standard_False;
464     TopExp_Explorer Explo(aFace, TopAbs_EDGE);
465     for (; Explo.More(); Explo.Next())
466     {
467       const TopoDS_Shape& anEdge = Explo.Current();
468       if (anEdge.IsSame(theEdge1))
469       {
470         e1found = Standard_True;
471         theOrOfE1OnFace = anEdge.Orientation();
472       }
473       if (anEdge.IsSame(theEdge2))
474       {
475         e2found = Standard_True;
476         theOrOfE2OnFace = anEdge.Orientation();
477       }
478       if (e1found && e2found)
479         return aFace;
480     }
481   }
482
483   TopoDS_Face NullFace;
484   return NullFace;
485 }
486
487 static Standard_Boolean FindClosestPoints(const TopoDS_Edge& theEdge1,
488                                           const TopoDS_Edge& theEdge2,
489                                           const TopTools_IndexedDataMapOfShapeListOfShape& theVFmap,
490                                           TopoDS_Face& theCommonFace,
491                                           Standard_Real& theMinSqDist,
492                                           TopAbs_Orientation& OrOfE1OnFace,
493                                           TopAbs_Orientation& OrOfE2OnFace,
494                                           Standard_Integer& theIndOnE1,
495                                           Standard_Integer& theIndOnE2,
496                                           gp_Pnt2d* PointsOnEdge1,
497                                           gp_Pnt2d* PointsOnEdge2)
498                                 
499 {
500   theCommonFace = FindCommonFace(theEdge1, theEdge2, theVFmap, OrOfE1OnFace, OrOfE2OnFace);
501   if (theCommonFace.IsNull())
502     return Standard_False;
503   
504   Standard_Real fpar1, lpar1, fpar2, lpar2;
505   Handle(Geom2d_Curve) PCurve1 = BRep_Tool::CurveOnSurface(theEdge1, theCommonFace, fpar1, lpar1);
506   Handle(Geom2d_Curve) PCurve2 = BRep_Tool::CurveOnSurface(theEdge2, theCommonFace, fpar2, lpar2);
507   PointsOnEdge1[0] = PCurve1->Value(fpar1);
508   PointsOnEdge1[1] = PCurve1->Value(lpar1);
509   PointsOnEdge2[0] = PCurve2->Value(fpar2);
510   PointsOnEdge2[1] = PCurve2->Value(lpar2);
511   theMinSqDist = RealLast();
512   theIndOnE1 = -1;
513   theIndOnE2 = -1;
514   for (Standard_Integer ind1 = 0; ind1 < 2; ind1++)
515     for (Standard_Integer ind2 = 0; ind2 < 2; ind2++)
516     {
517       Standard_Real aSqDist = PointsOnEdge1[ind1].SquareDistance(PointsOnEdge2[ind2]);
518       if (aSqDist < theMinSqDist)
519       {
520         theMinSqDist = aSqDist;
521         theIndOnE1 = ind1; theIndOnE2 = ind2;
522       }
523     }
524   return Standard_True;
525 }
526
527 //=======================================================================
528 //function : ReconstructMissedSeam
529 //purpose  : auxiliary
530 //=======================================================================
531 static void ReconstructMissedSeam(const TopTools_SequenceOfShape& theEdges,
532                                   const TopTools_SequenceOfShape& theRemovedEdges,
533                                   const TopTools_MapOfShape&      theUsedEdges,
534                                   const TopoDS_Face&              theFrefFace,
535                                   const TopoDS_Vertex&            theCurVertex,
536                                   const gp_Pnt2d&                 theCurPoint,
537                                   const Standard_Real             thePeriod,
538                                   const Standard_Real             theFaceCoordMin,
539                                   const Standard_Boolean          theIsU,
540                                   const Standard_Real             theCoordTol,
541                                   TopoDS_Edge&                    theNextEdge,
542                                   TopoDS_Wire&                    theNewWire,
543                                   gp_Pnt2d&                       theNextPoint,
544                                   gp_Pnt2d&                       theStartOfNextEdge,
545                                   TopoDS_Vertex&                  theLastVertexOfSeam,
546                                   TopTools_IndexedDataMapOfShapeListOfShape& theVEmap)
547                                   
548 {
549   Handle(Geom_Surface) RefSurf = BRep_Tool::Surface(theFrefFace);
550   GeomAbs_Shape aContinuity;
551   if (theIsU)
552     aContinuity = (RefSurf->IsUPeriodic())? GeomAbs_CN : GeomAbs_C0;
553   else
554     aContinuity = (RefSurf->IsVPeriodic())? GeomAbs_CN : GeomAbs_C0;
555
556   Standard_Integer IndCoord = theIsU? 1 : 2;
557   
558   Standard_Real SeamDir = 1.; //up or right
559   if (theIsU && Abs(theCurPoint.Coord(IndCoord) - theFaceCoordMin) <= theCoordTol)
560     SeamDir = -1.; //down
561   else if (!theIsU && Abs(theCurPoint.Coord(IndCoord) - theFaceCoordMin) > theCoordTol)
562     SeamDir = -1.; //left
563
564   //Consider as the candidate to be next edge:
565   //only the edge that has first point with X(or Y)-coordinate close to X(or Y)-coordinate of theCurPoint
566   //Choose from candidates the edge that is closest to theCurPoint in the defined direction SeamDir
567   Standard_Real MinDeltaSeamCoord = RealLast();
568   for (Standard_Integer ind = 1; ind <= theEdges.Length(); ind++)
569   {
570     const TopoDS_Edge& aCandidate = TopoDS::Edge(theEdges(ind));
571     if (theUsedEdges.Contains(aCandidate))
572       continue;
573     BRepAdaptor_Curve2d BAcurve2d(aCandidate, theFrefFace);
574     Standard_Real CandParam = (aCandidate.Orientation() == TopAbs_FORWARD)?
575       BAcurve2d.FirstParameter() : BAcurve2d.LastParameter();
576     gp_Pnt2d CandPoint = BAcurve2d.Value(CandParam);
577     Standard_Real DeltaCoord = Abs(CandPoint.Coord(IndCoord) - theCurPoint.Coord(IndCoord));
578     if (DeltaCoord > theCoordTol)
579       continue;
580     
581     Standard_Real DeltaSeamCoord = CandPoint.Coord(3-IndCoord) - theCurPoint.Coord(3-IndCoord);
582     DeltaSeamCoord *= SeamDir;
583     if (DeltaSeamCoord < 0.) //on the other side from CurPoint
584       continue;
585     if (DeltaSeamCoord < MinDeltaSeamCoord)
586     {
587       MinDeltaSeamCoord = DeltaSeamCoord;
588       theNextEdge = aCandidate;
589       theStartOfNextEdge = CandPoint;
590     }
591   }
592   
593   //Build missed seam edge
594   theLastVertexOfSeam = TopExp::FirstVertex(theNextEdge, Standard_True); //with orientation
595   TopoDS_Vertex V1, V2;
596   Standard_Real Param1, Param2, aCoord = 0.;
597   Handle(Geom_Curve) Iso;
598   
599   TopoDS_Edge aRemovedEdge; //try to find it in <RemovedEdges>
600   for (Standard_Integer i = 1; i <= theRemovedEdges.Length(); i++)
601   {
602     const TopoDS_Edge& anEdge = TopoDS::Edge(theRemovedEdges(i));
603     Handle(Geom2d_Curve) aPC = BRep_Tool::CurveOnSurface(anEdge, theFrefFace, Param1, Param2);
604     if (aPC.IsNull())
605       continue;
606     
607     GeomAbs_Shape aContOnRefFace = BRep_Tool::Continuity(anEdge, theFrefFace, theFrefFace);
608     if (aContOnRefFace > aContinuity)
609       aContinuity = aContOnRefFace;
610     
611     TopoDS_Vertex aV1, aV2;
612     TopExp::Vertices(anEdge, aV1, aV2);
613     if ((aV1.IsSame(theCurVertex) && aV2.IsSame(theLastVertexOfSeam)) ||
614         (aV1.IsSame(theLastVertexOfSeam) && aV2.IsSame(theCurVertex)))
615     {
616       aRemovedEdge = anEdge;
617       break;
618     }
619   }
620   if (aRemovedEdge.IsNull())
621   {
622     Standard_Real CurTol  = BRep_Tool::Tolerance(theCurVertex);
623     Standard_Real LastTol = BRep_Tool::Tolerance(theLastVertexOfSeam);
624     aCoord = (CurTol < LastTol)? theCurPoint.Coord(IndCoord) : theStartOfNextEdge.Coord(IndCoord);
625     Iso = (theIsU)? RefSurf->UIso(aCoord) : RefSurf->VIso(aCoord);
626     if (SeamDir > 0)
627     {
628       V1 = theCurVertex; V2 = theLastVertexOfSeam;
629       Param1 = theCurPoint.Coord(3-IndCoord); Param2 = theStartOfNextEdge.Coord(3-IndCoord);
630     }
631     else
632     {
633       V1 = theLastVertexOfSeam; V2 = theCurVertex;
634       Param1 = theStartOfNextEdge.Coord(3-IndCoord); Param2 = theCurPoint.Coord(3-IndCoord);
635     }
636   }
637   else
638   {
639     TopExp::Vertices(aRemovedEdge, V1, V2);
640     Iso = BRep_Tool::Curve(aRemovedEdge, Param1, Param2);
641   }
642   
643   TopoDS_Edge MissedSeam = BRepLib_MakeEdge(Iso, V1, V2, Param1, Param2);
644   BRep_Builder BB;
645   
646   //gp_Vec2d Offset(theUperiod, 0.);
647   gp_Vec2d Offset;
648   if (theIsU)
649     Offset.SetCoord(thePeriod, 0.);
650   else
651     Offset.SetCoord(0., thePeriod);
652   if (aRemovedEdge.IsNull())
653   {
654     Standard_Real SeamCoordOrigin = 0.;
655     //Correct Param1 and Param2 if needed:
656     //when iso-curve is periodic and Param1 and Param2 do not fit into SeamCoord-range of surface,
657     //(for example, V-range of sphere)
658     //BRepLib_MakeEdge may shift Param1 and Param2
659     Standard_Real InitialParam1 = Param1, InitialParam2 = Param2;
660     Handle(Geom_Curve) MissedCurve = BRep_Tool::Curve(MissedSeam, Param1, Param2);
661     if ((Param1 != InitialParam1 || Param2 != InitialParam2) &&
662         MissedCurve->IsPeriodic())
663     {
664       //Vorigin = -(MissedCurve->Period());
665       SeamCoordOrigin = -(Param1 - InitialParam1);
666     }
667     /////////////////////////////////////
668     //Handle(Geom2d_Line) PC1 = new Geom2d_Line(gp_Pnt2d(anU, Vorigin), gp_Dir2d(0., 1.));
669     Handle(Geom2d_Line) PC1;
670     if (theIsU)
671       PC1 = new Geom2d_Line(gp_Pnt2d(aCoord, SeamCoordOrigin), gp_Dir2d(0., 1.));
672     else
673       PC1 = new Geom2d_Line(gp_Pnt2d(SeamCoordOrigin, aCoord), gp_Dir2d(1., 0.));
674     
675     Handle(Geom2d_Curve) PC2 = Handle(Geom2d_Curve)::DownCast(PC1->Copy());
676     if (theIsU && SeamDir > 0)
677       Offset *= -1;
678     else if (!theIsU && SeamDir < 0)
679       Offset *= -1;
680     PC2->Translate(Offset);
681     
682     if (SeamDir > 0)
683       BB.UpdateEdge(MissedSeam, PC1, PC2, theFrefFace, 0.);
684     else
685       BB.UpdateEdge(MissedSeam, PC2, PC1, theFrefFace, 0.);
686
687     if (SeamDir < 0)
688       MissedSeam.Reverse();
689   }
690   else
691   {
692     TopoDS_Edge aSeam = aRemovedEdge;
693     aSeam.Orientation(TopAbs_FORWARD);
694     Handle(Geom2d_Curve) PC1 = BRep_Tool::CurveOnSurface(aSeam, theFrefFace, Param1, Param2);
695     aSeam.Reverse();
696     Handle(Geom2d_Curve) PC2 = BRep_Tool::CurveOnSurface(aSeam, theFrefFace, Param1, Param2);
697     Standard_Boolean IsSeam = (PC1 != PC2);
698     if (!IsSeam) //it was not a seam
699     {
700       aCoord = theCurPoint.Coord(IndCoord);
701       gp_Pnt2d PointOnRemovedEdge = PC1->Value(Param1);
702       Standard_Real CoordOfRemovededge = PointOnRemovedEdge.Coord(IndCoord);
703       if (Abs(aCoord - CoordOfRemovededge) > thePeriod/2)
704       {
705         Standard_Real Sign = (aCoord > CoordOfRemovededge)? 1 : -1;
706         Offset *= Sign;
707         PC1 = Handle(Geom2d_Curve)::DownCast(PC2->Copy());
708         PC1->Translate(Offset);
709       }
710       else
711       {
712         if (SeamDir > 0)
713           Offset *= -1;
714         PC2 = Handle(Geom2d_Curve)::DownCast(PC1->Copy());
715         PC2->Translate(Offset);
716       }
717     }
718     if (theCurVertex.IsSame(V1))
719       BB.UpdateEdge(MissedSeam, PC1, PC2, theFrefFace, 0.);
720     else
721     {
722       if (IsSeam)
723         BB.UpdateEdge(MissedSeam, PC1, PC2, theFrefFace, 0.);
724       else
725         BB.UpdateEdge(MissedSeam, PC2, PC1, theFrefFace, 0.);
726       
727       MissedSeam.Reverse();
728     }
729   }
730
731   BB.Continuity(MissedSeam, theFrefFace, theFrefFace, aContinuity);
732   BB.Add(theNewWire, MissedSeam);
733   //add newly created edge into VEmap
734   MissedSeam.Reverse();
735   TopExp::MapShapesAndUniqueAncestors(MissedSeam, TopAbs_VERTEX, TopAbs_EDGE, theVEmap);
736                   
737   BRepAdaptor_Curve2d BAcurve2d(theNextEdge, theFrefFace);
738   Standard_Real ParamOfNextPoint = (theNextEdge.Orientation() == TopAbs_FORWARD)?
739     BAcurve2d.LastParameter() : BAcurve2d.FirstParameter();
740   theNextPoint = BAcurve2d.Value(ParamOfNextPoint);
741 }
742
743 //=======================================================================
744 //function : SameSurf
745 //purpose  : auxiliary
746 //=======================================================================
747 static Standard_Boolean SameSurf(const Handle(Geom_Surface)& theS1, const Handle(Geom_Surface)& theS2)
748 {
749   static Standard_Real aCoefs[2] = { 0.3399811, 0.7745966 };
750
751   Standard_Real uf1, ul1, vf1, vl1, uf2, ul2, vf2, vl2;
752   theS1->Bounds(uf1, ul1, vf1, vl1);
753   theS2->Bounds(uf2, ul2, vf2, vl2);
754   Standard_Real aPTol = Precision::PConfusion();
755   if (Precision::IsNegativeInfinite(uf1))
756   {
757     if (!Precision::IsNegativeInfinite(uf2))
758     {
759       return Standard_False;
760     }
761     else
762     {
763       uf1 = Min(-1., (ul1 - 1.));
764     }
765   }
766   else
767   {
768     if (Precision::IsNegativeInfinite(uf2))
769     {
770       return Standard_False;
771     }
772     else
773     {
774       if (Abs(uf1 - uf2) > aPTol)
775       {
776         return Standard_False;
777       }
778     }
779   }
780   //
781   if (Precision::IsNegativeInfinite(vf1))
782   {
783     if (!Precision::IsNegativeInfinite(vf2))
784     {
785       return Standard_False;
786     }
787     else
788     {
789       vf1 = Min(-1., (vl1 - 1.));
790     }
791   }
792   else
793   {
794     if (Precision::IsNegativeInfinite(vf2))
795     {
796       return Standard_False;
797     }
798     else
799     {
800       if (Abs(vf1 - vf2) > aPTol)
801       {
802         return Standard_False;
803       }
804     }
805   }
806   //
807   if (Precision::IsPositiveInfinite(ul1))
808   {
809     if (!Precision::IsPositiveInfinite(ul2))
810     {
811       return Standard_False;
812     }
813     else
814     {
815       ul1 = Max(1., (uf1 + 1.));
816     }
817   }
818   else
819   {
820     if (Precision::IsPositiveInfinite(ul2))
821     {
822       return Standard_False;
823     }
824     else
825     {
826       if (Abs(ul1 - ul2) > aPTol)
827       {
828         return Standard_False;
829       }
830     }
831   }
832   //
833   if (Precision::IsPositiveInfinite(vl1))
834   {
835     if (!Precision::IsPositiveInfinite(vl2))
836     {
837       return Standard_False;
838     }
839     else
840     {
841       vl1 = Max(1., (vf1 + 1.));
842     }
843   }
844   else
845   {
846     if (Precision::IsPositiveInfinite(vl2))
847     {
848       return Standard_False;
849     }
850     else
851     {
852       if (Abs(vl1 - vl2) > aPTol)
853       {
854         return Standard_False;
855       }
856     }
857   }
858   //
859
860   Standard_Real u, v, du = (ul1 - uf1), dv = (vl1 - vf1);
861   Standard_Integer i, j;
862   for (i = 0; i < 2; ++i)
863   {
864     u = uf1 + aCoefs[i] * du;
865     for (j = 0; j < 2; ++j)
866     {
867       v = vf1 + aCoefs[j] * dv;
868       gp_Pnt aP1 = theS1->Value(u, v);
869       gp_Pnt aP2 = theS2->Value(u, v);
870       if (!aP1.IsEqual(aP2, aPTol))
871       {
872         return Standard_False;
873       }
874     }
875   }
876
877   return Standard_True;
878 }
879 //=======================================================================
880 //function : TransformPCurves
881 //purpose  : auxiliary
882 //=======================================================================
883 static void TransformPCurves(const TopoDS_Face& theRefFace,
884                              const TopoDS_Face& theFace,
885                              TopTools_MapOfShape& theMapEdgesWithTemporaryPCurves)
886 {
887   Handle(Geom_Surface) RefSurf = BRep_Tool::Surface(theRefFace);
888   if (RefSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
889     RefSurf = (Handle(Geom_RectangularTrimmedSurface)::DownCast(RefSurf))->BasisSurface();
890
891   Handle(Geom_Surface) SurfFace = BRep_Tool::Surface(theFace);
892   if (SurfFace->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
893     SurfFace = (Handle(Geom_RectangularTrimmedSurface)::DownCast(SurfFace))->BasisSurface();
894
895   Standard_Boolean ToModify = Standard_False,
896     ToTranslate = Standard_False,
897     ToRotate = Standard_False,
898     X_Reverse = Standard_False,
899     Y_Reverse = Standard_False,
900     ToProject = Standard_False;
901   
902   Standard_Real aTranslation = 0., anAngle = 0.;
903
904   //Get axes of surface of face and of surface of RefFace
905   Handle(Geom_ElementarySurface) ElemSurfFace = Handle(Geom_ElementarySurface)::DownCast(SurfFace);
906   Handle(Geom_ElementarySurface) ElemRefSurf = Handle(Geom_ElementarySurface)::DownCast(RefSurf);
907
908   if (!ElemSurfFace.IsNull() && !ElemRefSurf.IsNull())
909   {
910     gp_Ax3 AxisOfSurfFace = ElemSurfFace->Position();
911     gp_Ax3 AxisOfRefSurf = ElemRefSurf->Position();
912   
913     gp_Pnt OriginRefSurf  = AxisOfRefSurf.Location();
914
915     Standard_Real aParam = ElCLib::LineParameter(AxisOfSurfFace.Axis(), OriginRefSurf);
916
917     if (Abs(aParam) > Precision::PConfusion())
918       aTranslation = -aParam;
919
920     gp_Dir VdirSurfFace = AxisOfSurfFace.Direction();
921     gp_Dir VdirRefSurf  = AxisOfRefSurf.Direction();
922     gp_Dir XdirSurfFace = AxisOfSurfFace.XDirection();
923     gp_Dir XdirRefSurf  = AxisOfRefSurf.XDirection();
924   
925     gp_Dir CrossProd1 = AxisOfRefSurf.XDirection() ^ AxisOfRefSurf.YDirection();
926     gp_Dir CrossProd2 = AxisOfSurfFace.XDirection() ^ AxisOfSurfFace.YDirection();
927     if (CrossProd1 * CrossProd2 < 0.)
928       X_Reverse = Standard_True;
929
930     Standard_Real ScalProd = VdirSurfFace * VdirRefSurf;
931     if (ScalProd < 0.)
932       Y_Reverse = Standard_True;
933
934     if (!X_Reverse && !Y_Reverse)
935     {
936       gp_Dir DirRef = VdirRefSurf;
937       if (!AxisOfRefSurf.Direct())
938         DirRef.Reverse();
939       anAngle = XdirRefSurf.AngleWithRef(XdirSurfFace, DirRef);
940     }
941     else
942       anAngle = XdirRefSurf.Angle(XdirSurfFace);
943
944     ToRotate = (Abs(anAngle) > Precision::PConfusion());
945
946     ToTranslate = (Abs(aTranslation) > Precision::PConfusion());
947
948     ToModify = ToTranslate || ToRotate || X_Reverse || Y_Reverse;
949   }
950   else
951   {
952     if (!SameSurf(RefSurf, SurfFace))
953     {
954       ToProject = Standard_True;
955     }
956   }
957
958   BRep_Builder BB;
959   TopExp_Explorer Explo(theFace, TopAbs_EDGE);
960   for (; Explo.More(); Explo.Next())
961   {
962     const TopoDS_Edge& anEdge = TopoDS::Edge(Explo.Current());
963     if (BRep_Tool::Degenerated(anEdge) && ToModify)
964       continue;
965     if (BRepTools::IsReallyClosed(anEdge, theFace))
966       continue;
967
968     Standard_Real fpar, lpar;
969     Handle(Geom2d_Curve) PCurveOnRef = BRep_Tool::CurveOnSurface(anEdge, theRefFace, fpar, lpar);
970     if (!PCurveOnRef.IsNull() && !(ToModify || ToProject))
971       continue;
972
973     Handle(Geom2d_Curve) aPCurve = BRep_Tool::CurveOnSurface(anEdge, theFace, fpar, lpar);
974     Handle(Geom2d_Curve) aNewPCurve;
975     if (ToProject)
976     {
977       Handle(Geom_Curve) aC3d = BRep_Tool::Curve(anEdge, fpar, lpar);
978       aC3d = new Geom_TrimmedCurve(aC3d, fpar, lpar);
979       Standard_Real tol = BRep_Tool::Tolerance(anEdge);
980       tol = Min(tol, Precision::Approximation());
981       aNewPCurve =
982         GeomProjLib::Curve2d(aC3d, RefSurf);
983     }
984     else
985     {
986       aNewPCurve = Handle(Geom2d_Curve)::DownCast(aPCurve->Copy());
987     }
988     if (ToTranslate)
989       aNewPCurve->Translate(gp_Vec2d(0., aTranslation));
990     if (Y_Reverse)
991       aNewPCurve->Mirror(gp::OX2d());
992     if (X_Reverse)
993     {
994       aNewPCurve->Mirror(gp::OY2d());
995       aNewPCurve->Translate(gp_Vec2d(2*M_PI, 0.));
996     }
997     if (ToRotate)
998       aNewPCurve->Translate(gp_Vec2d(anAngle, 0.));
999
1000     theMapEdgesWithTemporaryPCurves.Add(anEdge);
1001     
1002     BB.UpdateEdge(anEdge, aNewPCurve, theRefFace, 0.);
1003     BB.Range(anEdge, fpar, lpar);
1004   }
1005 }
1006                              
1007 //=======================================================================
1008 //function : AddPCurves
1009 //purpose  : auxiliary
1010 //=======================================================================
1011 static void AddPCurves(const TopTools_SequenceOfShape& theFaces,
1012                        const TopoDS_Face&              theRefFace,
1013                        TopTools_MapOfShape& theMapEdgesWithTemporaryPCurves)
1014 {
1015   BRepAdaptor_Surface RefBAsurf(theRefFace, Standard_False);
1016
1017   GeomAbs_SurfaceType aType = RefBAsurf.GetType();
1018   if (aType == GeomAbs_Plane)
1019     return;
1020
1021   for (Standard_Integer i = 1; i <= theFaces.Length(); i++)
1022   {
1023     TopoDS_Face aFace = TopoDS::Face(theFaces(i));
1024     aFace.Orientation(TopAbs_FORWARD);
1025     if (aFace.IsSame(theRefFace))
1026       continue;
1027
1028     TransformPCurves(theRefFace, aFace, theMapEdgesWithTemporaryPCurves);
1029   }
1030 }
1031
1032 //=======================================================================
1033 //function : AddOrdinaryEdges
1034 //purpose  : auxiliary
1035 //=======================================================================
1036 // adds edges from the shape to the sequence
1037 // seams and equal edges are dropped
1038 // Returns true if one of original edges dropped
1039 static Standard_Boolean AddOrdinaryEdges(TopTools_SequenceOfShape& edges,
1040                                          const TopoDS_Shape aShape,
1041                                          Standard_Integer& anIndex,
1042                                          TopTools_SequenceOfShape& theRemovedEdges)
1043 {
1044   //map of edges
1045   TopTools_IndexedMapOfShape aNewEdges;
1046   //add edges without seams
1047   for(TopExp_Explorer exp(aShape,TopAbs_EDGE); exp.More(); exp.Next()) {
1048     TopoDS_Shape edge = exp.Current();
1049     if(aNewEdges.Contains(edge))
1050     {
1051       aNewEdges.RemoveKey(edge);
1052       theRemovedEdges.Append(edge);
1053     }
1054     else
1055       aNewEdges.Add(edge);
1056   }
1057
1058   Standard_Boolean isDropped = Standard_False;
1059   //merge edges and drop seams
1060   Standard_Integer i;
1061   for (i = 1; i <= edges.Length(); i++) {
1062     TopoDS_Shape current = edges(i);
1063     if(aNewEdges.Contains(current)) {
1064
1065       aNewEdges.RemoveKey(current);
1066       edges.Remove(i);
1067       theRemovedEdges.Append(current);
1068       i--;
1069
1070       if(!isDropped) {
1071         isDropped = Standard_True;
1072         anIndex = i;
1073       }
1074     }
1075   }
1076
1077   //add edges to the sequence
1078   for (i = 1; i <= aNewEdges.Extent(); i++)
1079     edges.Append(aNewEdges(i));
1080
1081   return isDropped;
1082 }
1083
1084 //=======================================================================
1085 //function : getCylinder
1086 //purpose  : auxiliary
1087 //=======================================================================
1088 static Standard_Boolean getCylinder(Handle(Geom_Surface)& theInSurface,
1089                                     gp_Cylinder& theOutCylinder)
1090 {
1091   Standard_Boolean isCylinder = Standard_False;
1092
1093   if (theInSurface->IsKind(STANDARD_TYPE(Geom_CylindricalSurface))) {
1094     Handle(Geom_CylindricalSurface) aGC = Handle(Geom_CylindricalSurface)::DownCast(theInSurface);
1095
1096     theOutCylinder = aGC->Cylinder();
1097     isCylinder = Standard_True;
1098   }
1099   else if (theInSurface->IsKind(STANDARD_TYPE(Geom_SurfaceOfRevolution))) {
1100     Handle(Geom_SurfaceOfRevolution) aRS =
1101       Handle(Geom_SurfaceOfRevolution)::DownCast(theInSurface);
1102     Handle(Geom_Curve) aBasis = aRS->BasisCurve();
1103     if (aBasis->IsKind(STANDARD_TYPE(Geom_Line))) {
1104       Handle(Geom_Line) aBasisLine = Handle(Geom_Line)::DownCast(aBasis);
1105       gp_Dir aDir = aRS->Direction();
1106       gp_Dir aBasisDir = aBasisLine->Position().Direction();
1107       if (aBasisDir.IsParallel(aDir, Precision::Angular())) {
1108         // basis line is parallel to the revolution axis: it is a cylinder
1109         gp_Pnt aLoc = aRS->Location();
1110         Standard_Real aR = aBasisLine->Lin().Distance(aLoc);
1111         gp_Ax3 aCylAx (aLoc, aDir);
1112
1113         theOutCylinder = gp_Cylinder(aCylAx, aR);
1114         isCylinder = Standard_True;
1115       }
1116     }
1117   }
1118   else if (theInSurface->IsKind(STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion))) {
1119     Handle(Geom_SurfaceOfLinearExtrusion) aLES =
1120       Handle(Geom_SurfaceOfLinearExtrusion)::DownCast(theInSurface);
1121     Handle(Geom_Curve) aBasis = aLES->BasisCurve();
1122     if (aBasis->IsKind(STANDARD_TYPE(Geom_Circle))) {
1123       Handle(Geom_Circle) aBasisCircle = Handle(Geom_Circle)::DownCast(aBasis);
1124       gp_Dir aDir = aLES->Direction();
1125       gp_Dir aBasisDir = aBasisCircle->Position().Direction();
1126       if (aBasisDir.IsParallel(aDir, Precision::Angular())) {
1127         // basis circle is normal to the extrusion axis: it is a cylinder
1128         gp_Pnt aLoc = aBasisCircle->Location();
1129         Standard_Real aR = aBasisCircle->Radius();
1130         gp_Ax3 aCylAx (aLoc, aDir);
1131
1132         theOutCylinder = gp_Cylinder(aCylAx, aR);
1133         isCylinder = Standard_True;
1134       }
1135     }
1136   }
1137   else {
1138   }
1139
1140   return isCylinder;
1141 }
1142
1143 //=======================================================================
1144 //function : ClearRts
1145 //purpose  : auxiliary
1146 //=======================================================================
1147 static Handle(Geom_Surface) ClearRts(const Handle(Geom_Surface)& aSurface)
1148 {
1149   if(aSurface->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
1150     Handle(Geom_RectangularTrimmedSurface) rts =
1151       Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurface);
1152     return rts->BasisSurface();
1153   }
1154   return aSurface;
1155 }
1156
1157 //=======================================================================
1158 //function : GetNormalToSurface
1159 //purpose  : Gets the normal to surface by the given parameter on edge.
1160 //           Returns True if normal was computed.
1161 //=======================================================================
1162 static Standard_Boolean GetNormalToSurface(const TopoDS_Face& theFace,
1163                                            const TopoDS_Edge& theEdge,
1164                                            const Standard_Real theP,
1165                                            gp_Dir& theNormal)
1166 {
1167   Standard_Real f, l;
1168   // get 2d curve to get point in 2d
1169   const Handle(Geom2d_Curve)& aC2d = BRep_Tool::CurveOnSurface(theEdge, theFace, f, l);
1170   if (aC2d.IsNull()) {
1171     return Standard_False;
1172   }
1173   //
1174   // 2d point
1175   gp_Pnt2d aP2d;
1176   aC2d->D0(theP, aP2d);
1177   //
1178   // get D1
1179   gp_Vec aDU, aDV;
1180   gp_Pnt aP3d;
1181   TopLoc_Location aLoc;
1182   const Handle(Geom_Surface)& aS = BRep_Tool::Surface(theFace, aLoc);
1183   aS->D1(aP2d.X(), aP2d.Y(), aP3d, aDU, aDV);
1184   //
1185   // compute normal
1186   gp_Vec aVNormal = aDU.Crossed(aDV);
1187   if (aVNormal.Magnitude() < Precision::Confusion()) {
1188     return Standard_False;
1189   }
1190   //
1191   if (theFace.Orientation() == TopAbs_REVERSED) {
1192     aVNormal.Reverse();
1193   }
1194   //
1195   aVNormal.Transform(aLoc.Transformation());
1196   theNormal = gp_Dir(aVNormal);
1197   return Standard_True;
1198 }
1199
1200 //=======================================================================
1201 //function : IsSameDomain
1202 //purpose  : 
1203 //=======================================================================
1204 static Standard_Boolean IsSameDomain(const TopoDS_Face& aFace,
1205                                      const TopoDS_Face& aCheckedFace,
1206                                      const Standard_Real theLinTol,
1207                                      const Standard_Real theAngTol,
1208                                      ShapeUpgrade_UnifySameDomain::DataMapOfFacePlane& theFacePlaneMap)
1209 {
1210   //checking the same handles
1211   TopLoc_Location L1, L2;
1212   Handle(Geom_Surface) S1, S2;
1213
1214   S1 = BRep_Tool::Surface(aFace,L1);
1215   S2 = BRep_Tool::Surface(aCheckedFace,L2);
1216
1217   if (S1 == S2 && L1 == L2)
1218     return Standard_True;
1219
1220   S1 = BRep_Tool::Surface(aFace);
1221   S2 = BRep_Tool::Surface(aCheckedFace);
1222
1223   S1 = ClearRts(S1);
1224   S2 = ClearRts(S2);
1225
1226   //Handle(Geom_OffsetSurface) aGOFS1, aGOFS2;
1227   //aGOFS1 = Handle(Geom_OffsetSurface)::DownCast(S1);
1228   //aGOFS2 = Handle(Geom_OffsetSurface)::DownCast(S2);
1229   //if (!aGOFS1.IsNull()) S1 = aGOFS1->BasisSurface();
1230   //if (!aGOFS2.IsNull()) S2 = aGOFS2->BasisSurface();
1231
1232   // case of two planar surfaces:
1233   // all kinds of surfaces checked, including b-spline and bezier
1234   GeomLib_IsPlanarSurface aPlanarityChecker1(S1, theLinTol);
1235   if (aPlanarityChecker1.IsPlanar()) {
1236     GeomLib_IsPlanarSurface aPlanarityChecker2(S2, theLinTol);
1237     if (aPlanarityChecker2.IsPlanar()) {
1238       gp_Pln aPln1 = aPlanarityChecker1.Plan();
1239       gp_Pln aPln2 = aPlanarityChecker2.Plan();
1240
1241       if (aPln1.Position().Direction().IsParallel(aPln2.Position().Direction(), theAngTol) &&
1242         aPln1.Distance(aPln2) < theLinTol)
1243       {
1244         Handle(Geom_Plane) aPlaneOfFaces;
1245         if (theFacePlaneMap.IsBound(aFace))
1246           aPlaneOfFaces = theFacePlaneMap(aFace);
1247         else if (theFacePlaneMap.IsBound(aCheckedFace))
1248           aPlaneOfFaces = theFacePlaneMap(aCheckedFace);
1249         else
1250           aPlaneOfFaces = new Geom_Plane(aPln1);
1251
1252         theFacePlaneMap.Bind(aFace, aPlaneOfFaces);
1253         theFacePlaneMap.Bind(aCheckedFace, aPlaneOfFaces);
1254         
1255         return Standard_True;
1256       }
1257     }
1258   }
1259
1260   // case of two elementary surfaces: use OCCT tool
1261   // elementary surfaces: ConicalSurface, CylindricalSurface,
1262   //                      Plane, SphericalSurface and ToroidalSurface
1263   if (S1->IsKind(STANDARD_TYPE(Geom_ElementarySurface)) &&
1264       S2->IsKind(STANDARD_TYPE(Geom_ElementarySurface)))
1265   {
1266     Handle(GeomAdaptor_Surface) aGA1 = new GeomAdaptor_Surface(S1);
1267     Handle(GeomAdaptor_Surface) aGA2 = new GeomAdaptor_Surface(S2);
1268
1269     Handle(BRepTopAdaptor_TopolTool) aTT1 = new BRepTopAdaptor_TopolTool();
1270     Handle(BRepTopAdaptor_TopolTool) aTT2 = new BRepTopAdaptor_TopolTool();
1271
1272     try {
1273       IntPatch_ImpImpIntersection anIIInt(aGA1, aTT1, aGA2, aTT2, theLinTol, theLinTol);
1274       if (!anIIInt.IsDone() || anIIInt.IsEmpty())
1275         return Standard_False;
1276
1277       return anIIInt.TangentFaces();
1278     }
1279     catch (Standard_Failure const&) {
1280       return Standard_False;
1281     }
1282   }
1283
1284   // case of two cylindrical surfaces, at least one of which is a swept surface
1285   // swept surfaces: SurfaceOfLinearExtrusion, SurfaceOfRevolution
1286   if ((S1->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)) ||
1287        S1->IsKind(STANDARD_TYPE(Geom_SweptSurface))) &&
1288       (S2->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)) ||
1289        S2->IsKind(STANDARD_TYPE(Geom_SweptSurface))))
1290   {
1291     gp_Cylinder aCyl1, aCyl2;
1292     if (getCylinder(S1, aCyl1) && getCylinder(S2, aCyl2)) {
1293       if (fabs(aCyl1.Radius() - aCyl2.Radius()) < theLinTol) {
1294         gp_Dir aDir1 = aCyl1.Position().Direction();
1295         gp_Dir aDir2 = aCyl2.Position().Direction();
1296         if (aDir1.IsParallel(aDir2, Precision::Angular())) {
1297           gp_Pnt aLoc1 = aCyl1.Location();
1298           gp_Pnt aLoc2 = aCyl2.Location();
1299           gp_Vec aVec12 (aLoc1, aLoc2);
1300           if (aVec12.SquareMagnitude() < theLinTol*theLinTol ||
1301               aVec12.IsParallel(aDir1, Precision::Angular())) {
1302             return Standard_True;
1303           }
1304         }
1305       }
1306     }
1307   }
1308
1309   return Standard_False;
1310 }
1311
1312 //=======================================================================
1313 //function : UpdateMapOfShapes
1314 //purpose  :
1315 //=======================================================================
1316 static void UpdateMapOfShapes(TopTools_MapOfShape& theMapOfShapes,
1317                               Handle(ShapeBuild_ReShape)& theContext)
1318 {
1319   for (TopTools_MapIteratorOfMapOfShape it(theMapOfShapes); it.More(); it.Next()) {
1320     const TopoDS_Shape& aShape = it.Value();
1321     TopoDS_Shape aContextShape = theContext->Apply(aShape);
1322     if (!aContextShape.IsSame(aShape))
1323       theMapOfShapes.Add(aContextShape);
1324   }
1325 }
1326
1327 //=======================================================================
1328 //function : GlueEdgesWithPCurves
1329 //purpose  : Glues the pcurves of the sequence of edges
1330 //           and glues their 3d curves
1331 //=======================================================================
1332 static TopoDS_Edge GlueEdgesWithPCurves(const TopTools_SequenceOfShape& aChain,
1333                                         const TopoDS_Vertex& FirstVertex,
1334                                         const TopoDS_Vertex& LastVertex)
1335 {
1336   Standard_Integer i, j;
1337
1338   TopoDS_Edge FirstEdge = TopoDS::Edge(aChain(1));
1339   TColGeom_SequenceOfSurface SurfSeq;
1340   NCollection_Sequence<TopLoc_Location> LocSeq;
1341   
1342   for (int aCurveIndex = 0;; aCurveIndex++)
1343   {
1344     Handle(Geom2d_Curve) aCurve;
1345     Handle(Geom_Surface) aSurface;
1346     TopLoc_Location aLocation;
1347     Standard_Real aFirst, aLast;
1348     BRep_Tool::CurveOnSurface (FirstEdge, aCurve, aSurface, aLocation, aFirst, aLast, aCurveIndex);
1349     if (aCurve.IsNull())
1350       break;
1351
1352     SurfSeq.Append(aSurface);
1353     LocSeq.Append(aLocation);
1354   }
1355
1356   Standard_Real fpar, lpar;
1357   BRep_Tool::Range(FirstEdge, fpar, lpar);
1358   TopoDS_Edge PrevEdge = FirstEdge;
1359   TopoDS_Vertex CV;
1360   Standard_Real MaxTol = 0.;
1361   
1362   TopoDS_Edge ResEdge;
1363   BRep_Builder BB;
1364
1365   Standard_Integer nb_curve = aChain.Length();   //number of curves
1366   TColGeom_Array1OfBSplineCurve tab_c3d(0,nb_curve-1);                    //array of the curves
1367   TColStd_Array1OfReal tabtolvertex(0,nb_curve-1); //(0,nb_curve-2);  //array of the tolerances
1368     
1369   TopoDS_Vertex PrevVertex = FirstVertex;
1370   for (i = 1; i <= nb_curve; i++)
1371   {
1372     TopoDS_Edge anEdge = TopoDS::Edge(aChain(i));
1373     TopoDS_Vertex VF, VL;
1374     TopExp::Vertices(anEdge, VF, VL);
1375     Standard_Boolean ToReverse = (!VF.IsSame(PrevVertex));
1376     
1377     Standard_Real Tol1 = BRep_Tool::Tolerance(VF);
1378     Standard_Real Tol2 = BRep_Tool::Tolerance(VL);
1379     if (Tol1 > MaxTol)
1380       MaxTol = Tol1;
1381     if (Tol2 > MaxTol)
1382       MaxTol = Tol2;
1383     
1384     if (i > 1)
1385     {
1386       TopExp::CommonVertex(PrevEdge, anEdge, CV);
1387       Standard_Real Tol = BRep_Tool::Tolerance(CV);
1388       tabtolvertex(i-2) = Tol;
1389     }
1390     
1391     Handle(Geom_Curve) aCurve = BRep_Tool::Curve(anEdge, fpar, lpar);
1392     Handle(Geom_TrimmedCurve) aTrCurve = new Geom_TrimmedCurve(aCurve, fpar, lpar);
1393     tab_c3d(i-1) = GeomConvert::CurveToBSplineCurve(aTrCurve);
1394     GeomConvert::C0BSplineToC1BSplineCurve(tab_c3d(i-1), Precision::Confusion());
1395     if (ToReverse)
1396       tab_c3d(i-1)->Reverse();
1397     PrevVertex = (ToReverse)? VF : VL;
1398     PrevEdge = anEdge;
1399   }
1400   Handle(TColGeom_HArray1OfBSplineCurve)  concatcurve;     //array of the concatenated curves
1401   Handle(TColStd_HArray1OfInteger)        ArrayOfIndices;  //array of the remaining Vertex
1402   Standard_Boolean closed_flag = Standard_False;
1403   GeomConvert::ConcatC1(tab_c3d,
1404                         tabtolvertex,
1405                         ArrayOfIndices,
1406                         concatcurve,
1407                         closed_flag,
1408                         Precision::Confusion());   //C1 concatenation
1409   
1410   if (concatcurve->Length() > 1)
1411   {
1412     GeomConvert_CompCurveToBSplineCurve Concat(concatcurve->Value(concatcurve->Lower()));
1413     
1414     for (i = concatcurve->Lower()+1; i <= concatcurve->Upper(); i++)
1415       Concat.Add( concatcurve->Value(i), MaxTol, Standard_True );
1416     
1417     concatcurve->SetValue(concatcurve->Lower(), Concat.BSplineCurve());
1418   }
1419   Handle(Geom_BSplineCurve) ResCurve = concatcurve->Value(concatcurve->Lower());
1420   
1421   TColGeom2d_SequenceOfBoundedCurve ResPCurves;
1422   for (j = 1; j <= SurfSeq.Length(); j++)
1423   {
1424     TColGeom2d_Array1OfBSplineCurve tab_c2d(0,nb_curve-1); //array of the pcurves
1425     
1426     PrevVertex = FirstVertex;
1427     PrevEdge = FirstEdge;
1428     for (i = 1; i <= nb_curve; i++)
1429     {
1430       TopoDS_Edge anEdge = TopoDS::Edge(aChain(i));
1431       TopoDS_Vertex VF, VL;
1432       TopExp::Vertices(anEdge, VF, VL);
1433       Standard_Boolean ToReverse = (!VF.IsSame(PrevVertex));
1434
1435       Handle(Geom2d_Curve) aPCurve =
1436         BRep_Tool::CurveOnSurface(anEdge, SurfSeq(j), LocSeq(j), fpar, lpar);
1437       if (aPCurve.IsNull())
1438         continue;
1439       Handle(Geom2d_TrimmedCurve) aTrPCurve = new Geom2d_TrimmedCurve(aPCurve, fpar, lpar);
1440       tab_c2d(i-1) = Geom2dConvert::CurveToBSplineCurve(aTrPCurve);
1441       Geom2dConvert::C0BSplineToC1BSplineCurve(tab_c2d(i-1), Precision::Confusion());
1442       if (ToReverse)
1443         tab_c2d(i-1)->Reverse();
1444       PrevVertex = (ToReverse)? VF : VL;
1445       PrevEdge = anEdge;
1446     }
1447     Handle(TColGeom2d_HArray1OfBSplineCurve)  concatc2d;     //array of the concatenated curves
1448     Handle(TColStd_HArray1OfInteger)        ArrayOfInd2d;  //array of the remaining Vertex
1449     closed_flag = Standard_False;
1450     Geom2dConvert::ConcatC1(tab_c2d,
1451                             tabtolvertex,
1452                             ArrayOfInd2d,
1453                             concatc2d,
1454                             closed_flag,
1455                             Precision::Confusion());   //C1 concatenation
1456     
1457     if (concatc2d->Length() > 1)
1458     {
1459       Geom2dConvert_CompCurveToBSplineCurve Concat2d(concatc2d->Value(concatc2d->Lower()));
1460       
1461       for (i = concatc2d->Lower()+1; i <= concatc2d->Upper(); i++)
1462         Concat2d.Add( concatc2d->Value(i), MaxTol, Standard_True );
1463       
1464       concatc2d->SetValue(concatc2d->Lower(), Concat2d.BSplineCurve());
1465     }
1466     Handle(Geom2d_BSplineCurve) aResPCurve = concatc2d->Value(concatc2d->Lower());
1467     ResPCurves.Append(aResPCurve);
1468   }
1469   
1470   ResEdge = BRepLib_MakeEdge(ResCurve,
1471                              FirstVertex, LastVertex,
1472                              ResCurve->FirstParameter(), ResCurve->LastParameter());
1473   BB.SameRange(ResEdge, Standard_False);
1474   BB.SameParameter(ResEdge, Standard_False);
1475   for (j = 1; j <= ResPCurves.Length(); j++)
1476   {
1477     BB.UpdateEdge(ResEdge, ResPCurves(j), SurfSeq(j), LocSeq(j), MaxTol);
1478     BB.Range(ResEdge, SurfSeq(j), LocSeq(j), ResPCurves(j)->FirstParameter(), ResPCurves(j)->LastParameter());
1479   }
1480
1481   BRepLib::SameParameter(ResEdge, MaxTol, Standard_True);
1482   
1483   return ResEdge;
1484 }
1485
1486 //=======================================================================
1487 //function : UnionPCurves
1488 //purpose  : 
1489 //=======================================================================
1490
1491 void ShapeUpgrade_UnifySameDomain::UnionPCurves(const TopTools_SequenceOfShape& theChain,
1492                                                 TopoDS_Edge& theEdge)
1493 {
1494   Standard_Real aFirst3d, aLast3d;
1495   Handle(Geom_Curve) aCurve = BRep_Tool::Curve (theEdge, aFirst3d, aLast3d);
1496   Standard_Real aTolEdge = BRep_Tool::Tolerance(theEdge);
1497   Standard_Real aMaxTol = aTolEdge;
1498   
1499   TopTools_SequenceOfShape aFaceSeq;
1500
1501   const TopoDS_Edge& aFirstEdge = TopoDS::Edge(theChain.Value(1));
1502   const TopTools_ListOfShape& aFaceList = myEFmap.FindFromKey (aFirstEdge);
1503   TopTools_ListIteratorOfListOfShape anItl (aFaceList);
1504   for (; anItl.More(); anItl.Next())
1505   {
1506     TopoDS_Face aFace = TopoDS::Face (anItl.Value());
1507     if (myFacePlaneMap.IsBound(aFace))
1508       continue;
1509
1510     if (myFaceNewFace.IsBound(aFace))
1511       aFace = TopoDS::Face (myFaceNewFace(aFace));
1512
1513     BRepAdaptor_Surface aBAsurf (aFace, Standard_False);
1514     if (aBAsurf.GetType() == GeomAbs_Plane)
1515       continue;
1516
1517     TopLoc_Location aLoc;
1518     Handle(Geom_Surface) aSurf = BRep_Tool::Surface (aFace, aLoc);
1519     Standard_Boolean isFound = Standard_False;
1520     for (Standard_Integer ii = 1; ii <= aFaceSeq.Length(); ii++)
1521     {
1522       TopLoc_Location aPrevLoc;
1523       Handle(Geom_Surface) aPrevSurf = BRep_Tool::Surface (TopoDS::Face(aFaceSeq(ii)), aPrevLoc);
1524       if (aPrevSurf == aSurf && aPrevLoc == aLoc)
1525       {
1526         isFound = Standard_True;
1527         break;
1528       }
1529     }
1530     if (isFound)
1531       continue;
1532     
1533     Standard_Real aFirst, aLast;
1534     Handle(Geom2d_Curve) aPCurve = BRep_Tool::CurveOnSurface (aFirstEdge, aFace, aFirst, aLast);
1535
1536     aFaceSeq.Append (aFace);
1537   }
1538   
1539   TColGeom2d_SequenceOfCurve ResPCurves;
1540   TColStd_SequenceOfReal ResFirsts;
1541   TColStd_SequenceOfReal ResLasts;
1542   TColStd_SequenceOfReal aTolVerSeq;
1543   TopoDS_Edge aPrevEdge;
1544
1545   for (Standard_Integer j = 1; j <= aFaceSeq.Length(); j++)
1546   {
1547     TColGeom2d_SequenceOfCurve aPCurveSeq;
1548     TColStd_SequenceOfReal aFirstsSeq;
1549     TColStd_SequenceOfReal aLastsSeq;
1550     TColStd_SequenceOfBoolean aForwardsSeq;
1551     GeomAbs_CurveType aCurrentType = GeomAbs_OtherCurve;
1552
1553     Standard_Real aFirst, aLast;
1554     for (Standard_Integer i = 1; i <= theChain.Length(); i++)
1555     {
1556       TopoDS_Edge anEdge = TopoDS::Edge(theChain.Value(i));
1557       Standard_Boolean isForward = (anEdge.Orientation() != TopAbs_REVERSED);
1558
1559       Handle(Geom2d_Curve) aPCurve = BRep_Tool::CurveOnSurface(anEdge, TopoDS::Face(aFaceSeq(j)),
1560                                                                aFirst, aLast);
1561       if (aPCurve.IsNull())
1562         continue;
1563
1564       Geom2dAdaptor_Curve anAdaptor(aPCurve);
1565       GeomAbs_CurveType aType = anAdaptor.GetType();
1566
1567       Handle(Geom2d_Line) aLine;
1568       if (aType == GeomAbs_BSplineCurve ||
1569           aType == GeomAbs_BezierCurve)
1570         TryMakeLine (aPCurve, aFirst, aLast, aLine);
1571       if (!aLine.IsNull())
1572       {
1573         aPCurve = aLine;
1574         anAdaptor.Load (aPCurve);
1575         aType = GeomAbs_Line;
1576       }
1577
1578       if (aPCurveSeq.IsEmpty()) {
1579         Handle(Geom2d_Curve) aCopyPCurve = Handle(Geom2d_Curve)::DownCast(aPCurve->Copy());
1580         aPCurveSeq.Append(aCopyPCurve);
1581         aFirstsSeq.Append(aFirst);
1582         aLastsSeq.Append(aLast);
1583         aForwardsSeq.Append(isForward);
1584         aCurrentType = aType;
1585         aPrevEdge = anEdge;
1586         continue;
1587       }
1588       
1589       Standard_Boolean isSameCurve = Standard_False;
1590       Standard_Real aNewF = aFirst;
1591       Standard_Real aNewL = aLast;
1592       if (aPCurve == aPCurveSeq.Last())
1593       {
1594         isSameCurve = Standard_True;
1595       }
1596       else if (aType == aCurrentType)
1597       {
1598         Geom2dAdaptor_Curve aPrevAdaptor(aPCurveSeq.Last());
1599         switch (aType) {
1600         case GeomAbs_Line: {
1601           gp_Lin2d aPrevLin = aPrevAdaptor.Line();
1602           gp_Pnt2d aFirstP2d = aPCurve->Value (aFirst);
1603           gp_Pnt2d aLastP2d  = aPCurve->Value (aLast);
1604           if (aPrevLin.Contains (aFirstP2d,  Precision::Confusion()) &&
1605               aPrevLin.Contains (aLastP2d,  Precision::Confusion()))
1606           {
1607             isSameCurve = Standard_True;
1608             gp_Pnt2d p1 = anAdaptor.Value(aFirst);
1609             gp_Pnt2d p2 = anAdaptor.Value(aLast);
1610             aNewF = ElCLib::Parameter(aPrevLin, p1);
1611             aNewL = ElCLib::Parameter(aPrevLin, p2);
1612             if (aNewF > aNewL)
1613             {
1614               Standard_Real aTmp = aNewF;
1615               aNewF = aNewL;
1616               aNewL = aTmp;
1617             }
1618           }
1619           break;
1620         }
1621         case GeomAbs_Circle: {
1622           gp_Circ2d aCirc = anAdaptor.Circle();
1623           gp_Circ2d aPrevCirc = aPrevAdaptor.Circle();
1624           if (aCirc.Location().Distance(aPrevCirc.Location()) <= Precision::Confusion() &&
1625               Abs(aCirc.Radius() - aPrevCirc.Radius()) <= Precision::Confusion())
1626           {
1627             isSameCurve = Standard_True;
1628             gp_Pnt2d p1 = anAdaptor.Value(aFirst);
1629             gp_Pnt2d p2 = anAdaptor.Value(aLast);
1630             aNewF = ElCLib::Parameter(aPrevCirc, p1);
1631             aNewL = ElCLib::Parameter(aPrevCirc, p2);
1632             if (aNewF > aNewL)
1633             {
1634               Standard_Real aTmp = aNewF;
1635               aNewF = aNewL;
1636               aNewL = aTmp;
1637             }
1638           }
1639           break;
1640         }
1641         default:
1642           break;
1643         }
1644       }
1645       if (isSameCurve) {
1646         if (aForwardsSeq.Last() == Standard_True)
1647           aLastsSeq.ChangeLast() = aNewL;
1648         else
1649           aFirstsSeq.ChangeLast() = aNewF;
1650       }
1651       else
1652       {
1653         Handle(Geom2d_Curve) aCopyPCurve = Handle(Geom2d_Curve)::DownCast(aPCurve->Copy());
1654         aPCurveSeq.Append(aCopyPCurve);
1655         aFirstsSeq.Append(aFirst);
1656         aLastsSeq.Append(aLast);
1657         aForwardsSeq.Append(isForward);
1658         aCurrentType = aType;
1659         TopoDS_Vertex aV;
1660         TopExp::CommonVertex(aPrevEdge, anEdge, aV);
1661         Standard_Real aTol = BRep_Tool::Tolerance(aV);
1662         aTolVerSeq.Append (aTol);
1663       }
1664       aPrevEdge = anEdge;
1665     }
1666
1667     Handle(Geom2d_Curve) aResPCurve;
1668     Standard_Real aResFirst, aResLast;
1669     if (aPCurveSeq.Length() == 1) {
1670       aResPCurve = aPCurveSeq.Last();
1671       aResFirst = aFirstsSeq.Last();
1672       aResLast = aLastsSeq.Last();
1673       if (aForwardsSeq.Last() == Standard_False)
1674       {
1675         Standard_Real aNewLast  = aResPCurve->ReversedParameter (aResFirst);
1676         Standard_Real aNewFirst = aResPCurve->ReversedParameter (aResLast);
1677         aResPCurve->Reverse();
1678         aResFirst = aNewFirst;
1679         aResLast  = aNewLast;
1680       }
1681     }
1682     else {
1683       //C1 concatenation for PCurveSeq
1684       TColGeom2d_Array1OfBSplineCurve tab_c2d(0, aPCurveSeq.Length() - 1);
1685       for (Standard_Integer i = 1; i <= aPCurveSeq.Length(); i++) {
1686         Handle(Geom2d_TrimmedCurve) aTrPCurve = new Geom2d_TrimmedCurve(aPCurveSeq(i), aFirstsSeq(i), aLastsSeq(i));
1687         if (aForwardsSeq(i) == Standard_False)
1688         {
1689           aTrPCurve->Reverse();
1690         }
1691         tab_c2d(i - 1) = Geom2dConvert::CurveToBSplineCurve(aTrPCurve);
1692         Geom2dConvert::C0BSplineToC1BSplineCurve(tab_c2d(i - 1), Precision::Confusion());
1693       }
1694
1695       TColStd_Array1OfReal tabtolvertex(0, aTolVerSeq.Length() - 1);
1696       for (Standard_Integer i = 1; i <= aTolVerSeq.Length(); i++)
1697       {
1698         Standard_Real aTol = aTolVerSeq(i);
1699         tabtolvertex(i - 1) = aTol;
1700         if (aTol > aMaxTol)
1701           aMaxTol = aTol;
1702       }
1703
1704       Handle(TColGeom2d_HArray1OfBSplineCurve)  concatc2d;     //array of the concatenated curves
1705       Handle(TColStd_HArray1OfInteger)        ArrayOfInd2d;  //array of the remaining Vertex
1706       Standard_Boolean aClosedFlag = Standard_False;
1707       Geom2dConvert::ConcatC1(tab_c2d,
1708         tabtolvertex,
1709         ArrayOfInd2d,
1710         concatc2d,
1711         aClosedFlag,
1712         Precision::Confusion());   //C1 concatenation
1713
1714       if (concatc2d->Length() > 1)
1715       {
1716         Geom2dConvert_CompCurveToBSplineCurve Concat2d(concatc2d->Value(concatc2d->Lower()));
1717
1718         for (Standard_Integer i = concatc2d->Lower() + 1; i <= concatc2d->Upper(); i++)
1719           Concat2d.Add(concatc2d->Value(i), aMaxTol, Standard_True);
1720
1721         concatc2d->SetValue(concatc2d->Lower(), Concat2d.BSplineCurve());
1722       }
1723       Handle(Geom2d_BSplineCurve) aBSplineCurve = concatc2d->Value(concatc2d->Lower());
1724       aResPCurve = aBSplineCurve;
1725       aResFirst = aBSplineCurve->FirstParameter();
1726       aResLast = aBSplineCurve->LastParameter();
1727     }
1728     ResPCurves.Append(aResPCurve);
1729     ResFirsts.Append(aResFirst);
1730     ResLasts.Append(aResLast);
1731   }
1732
1733   BRep_Builder aBuilder;
1734
1735   //Check the results for consistency
1736   Standard_Boolean IsBadRange = Standard_False;
1737   Standard_Real aRange3d = aLast3d - aFirst3d;
1738   for (Standard_Integer ii = 1; ii <= ResPCurves.Length(); ii++)
1739   {
1740     Standard_Real aRange = ResLasts(ii) - ResFirsts(ii);
1741     if (Abs (aRange3d - aRange) > aMaxTol)
1742       IsBadRange = Standard_True;
1743   }
1744   
1745   if (IsBadRange)
1746   {
1747     for (Standard_Integer ii = 1; ii <= ResPCurves.Length(); ii++)
1748     {
1749       const TopoDS_Face& aFace = TopoDS::Face (aFaceSeq(ii));
1750       Handle(Geom_Surface) aSurf = BRep_Tool::Surface (aFace);
1751       Handle(ShapeAnalysis_Surface) aSAS = new ShapeAnalysis_Surface (aSurf);
1752       ShapeConstruct_ProjectCurveOnSurface aToolProj;
1753       aToolProj.Init (aSAS, Precision::Confusion());
1754       Handle(Geom2d_Curve) aNewPCurve;
1755       if (aToolProj.Perform(aCurve, aFirst3d, aLast3d, aNewPCurve))
1756         ResPCurves(ii) = aNewPCurve;
1757       else
1758       {
1759         //Reparametrize pcurve
1760         Handle(Geom2d_TrimmedCurve) aTrPCurve =
1761           new Geom2d_TrimmedCurve (ResPCurves(ii), ResFirsts(ii),  ResLasts(ii));
1762         Handle(Geom2d_BSplineCurve) aBSplinePCurve = Geom2dConvert::CurveToBSplineCurve(aTrPCurve);
1763         TColStd_Array1OfReal aKnots (1, aBSplinePCurve->NbKnots());
1764         aBSplinePCurve->Knots (aKnots);
1765         BSplCLib::Reparametrize (aFirst3d, aLast3d, aKnots);
1766         aBSplinePCurve->SetKnots (aKnots);
1767         ResPCurves(ii) = aBSplinePCurve;
1768       }
1769       ResFirsts(ii) = aFirst3d;
1770       ResLasts(ii)  = aLast3d;
1771     }
1772   }
1773
1774   //Reparametrize 3d curve if needed
1775   if (!ResPCurves.IsEmpty())
1776   {
1777     if (Abs (aFirst3d - ResFirsts(1)) > aMaxTol ||
1778         Abs (aLast3d  - ResLasts(1))  > aMaxTol)
1779     {
1780       GeomAdaptor_Curve aGAcurve (aCurve);
1781       GeomAbs_CurveType aType = aGAcurve.GetType();
1782       if (aType == GeomAbs_Line)
1783       {
1784         gp_Lin aLin = aGAcurve.Line();
1785         gp_Dir aDir = aLin.Direction();
1786         gp_Pnt aPnt = aGAcurve.Value (aFirst3d);
1787         gp_Vec anOffset = -aDir;
1788         anOffset *= ResFirsts(1);
1789         aPnt.Translate (anOffset);
1790         Handle(Geom_Line) aLine = new Geom_Line (aPnt, aDir);
1791         aBuilder.UpdateEdge (theEdge, aLine, aTolEdge);
1792         aBuilder.Range(theEdge, ResFirsts(1), ResLasts(1));
1793       }
1794       else if (aType == GeomAbs_Circle)
1795       {
1796         gp_Circ aCirc = aGAcurve.Circle();
1797         Standard_Real aRadius = aCirc.Radius();
1798         gp_Ax2 aPosition = aCirc.Position();
1799         gp_Ax1 anAxis = aPosition.Axis();
1800         Standard_Real anOffset = aFirst3d - ResFirsts(1);
1801         aPosition.Rotate (anAxis, anOffset);
1802         Handle(Geom_Circle) aCircle = new Geom_Circle (aPosition, aRadius);
1803         aBuilder.UpdateEdge (theEdge, aCircle, aTolEdge);
1804         aBuilder.Range(theEdge, ResFirsts(1), ResLasts(1));
1805       }
1806       else //general case
1807       {
1808         for (Standard_Integer ii = 1; ii <= ResPCurves.Length(); ii++)
1809         {
1810           if (Abs (aFirst3d - ResFirsts(ii)) > Precision::Confusion() ||
1811               Abs (aLast3d  - ResLasts(ii))  > Precision::Confusion())
1812           {
1813             Handle(Geom2d_TrimmedCurve) aTrPCurve =
1814               new Geom2d_TrimmedCurve (ResPCurves(ii), ResFirsts(ii),  ResLasts(ii));
1815             Handle(Geom2d_BSplineCurve) aBSplinePCurve = Geom2dConvert::CurveToBSplineCurve(aTrPCurve);
1816             TColStd_Array1OfReal aKnots (1, aBSplinePCurve->NbKnots());
1817             aBSplinePCurve->Knots (aKnots);
1818             BSplCLib::Reparametrize (aFirst3d, aLast3d, aKnots);
1819             aBSplinePCurve->SetKnots (aKnots);
1820             ResPCurves(ii) = aBSplinePCurve;
1821           }
1822         }
1823       }
1824     }
1825   }
1826
1827   for (Standard_Integer j = 1; j <= ResPCurves.Length(); j++)
1828   {
1829     aBuilder.UpdateEdge(theEdge, ResPCurves(j), TopoDS::Face(aFaceSeq(j)), aTolEdge);
1830   }
1831 }
1832
1833 //=======================================================================
1834 //function : MergeSubSeq
1835 //purpose  : Merges a sequence of edges into one edge if possible
1836 //=======================================================================
1837
1838 Standard_Boolean ShapeUpgrade_UnifySameDomain::MergeSubSeq(const TopTools_SequenceOfShape& theChain,
1839                                                            const TopTools_IndexedDataMapOfShapeListOfShape& theVFmap,
1840                                                            TopoDS_Edge& OutEdge)
1841 {
1842   ShapeAnalysis_Edge sae;
1843   BRep_Builder B;
1844   // union edges in chain
1845   int j;
1846   Standard_Real fp1,lp1,fp2,lp2;
1847   Standard_Boolean IsUnionOfLinesPossible = Standard_True;
1848   Standard_Boolean IsUnionOfCirclesPossible = Standard_True;
1849   Handle(Geom_Curve) c3d1, c3d2;
1850   for(j = 1; j < theChain.Length(); j++) 
1851   {
1852     TopoDS_Edge edge1 = TopoDS::Edge(theChain.Value(j));
1853     TopoDS_Edge edge2 = TopoDS::Edge(theChain.Value(j+1));
1854
1855     if (BRep_Tool::Degenerated(edge1) &&
1856         BRep_Tool::Degenerated(edge2))
1857     {
1858       //Find the closest points in 2d
1859       TopoDS_Edge edgeFirst = TopoDS::Edge(theChain.First());
1860       TopoDS_Edge edgeLast  = TopoDS::Edge(theChain.Last());
1861       TopoDS_Face CommonFace;
1862       Standard_Real MinSqDist;
1863       TopAbs_Orientation OrOfE1OnFace, OrOfE2OnFace;
1864       Standard_Integer IndOnE1, IndOnE2;
1865       gp_Pnt2d PointsOnEdge1 [2], PointsOnEdge2 [2];
1866       if (!FindClosestPoints(edgeFirst, edgeLast, theVFmap, CommonFace,
1867                              MinSqDist, OrOfE1OnFace, OrOfE2OnFace,
1868                              IndOnE1, IndOnE2, PointsOnEdge1, PointsOnEdge2))
1869         return Standard_False;
1870       
1871       //Define indices corresponding to extremities of future edge
1872       IndOnE1 = 1 - IndOnE1;
1873       IndOnE2 = 1 - IndOnE2;
1874
1875       //Construct new degenerated edge
1876       gp_Pnt2d StartPoint = PointsOnEdge1[IndOnE1];
1877       gp_Pnt2d EndPoint   = PointsOnEdge2[IndOnE2];
1878       if ((OrOfE1OnFace == TopAbs_FORWARD  && IndOnE1 == 1) ||
1879           (OrOfE1OnFace == TopAbs_REVERSED && IndOnE1 == 0))
1880       { gp_Pnt2d Tmp = StartPoint; StartPoint = EndPoint; EndPoint = Tmp; }
1881       
1882       Handle(Geom2d_Line) aLine = GCE2d_MakeLine(StartPoint, EndPoint);
1883
1884       TopoDS_Vertex aVertex = TopExp::FirstVertex(edgeFirst);
1885       TopoDS_Vertex StartVertex = aVertex, EndVertex = aVertex;
1886       StartVertex.Orientation(TopAbs_FORWARD);
1887       EndVertex.Orientation(TopAbs_REVERSED);
1888       
1889       TopoDS_Edge NewEdge;
1890       B.MakeEdge(NewEdge);
1891       B.UpdateEdge(NewEdge, aLine, CommonFace, Precision::Confusion());
1892       B.Range(NewEdge, 0., StartPoint.Distance(EndPoint));
1893       B.Add (NewEdge, StartVertex);
1894       B.Add (NewEdge, EndVertex);
1895       B.Degenerated(NewEdge, Standard_True);
1896       OutEdge = NewEdge;
1897       return Standard_True;
1898     }
1899     
1900     c3d1 = BRep_Tool::Curve(edge1,fp1,lp1);
1901     c3d2 = BRep_Tool::Curve(edge2,fp2,lp2);
1902
1903     if(c3d1.IsNull() || c3d2.IsNull()) 
1904       return Standard_False;
1905
1906     if (c3d1->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
1907       Handle(Geom_TrimmedCurve) tc =
1908         Handle(Geom_TrimmedCurve)::DownCast(c3d1);
1909       c3d1 = tc->BasisCurve();
1910     }
1911     if (c3d2->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
1912       Handle(Geom_TrimmedCurve) tc =
1913         Handle(Geom_TrimmedCurve)::DownCast(c3d2);
1914       c3d2 = tc->BasisCurve();
1915     }
1916     if( c3d1->IsKind(STANDARD_TYPE(Geom_Line)) && c3d2->IsKind(STANDARD_TYPE(Geom_Line)) ) {
1917       Handle(Geom_Line) L1 = Handle(Geom_Line)::DownCast(c3d1);
1918       Handle(Geom_Line) L2 = Handle(Geom_Line)::DownCast(c3d2);
1919       gp_Dir Dir1 = L1->Position().Direction();
1920       gp_Dir Dir2 = L2->Position().Direction();
1921       if (!Dir1.IsParallel (Dir2, myAngTol))  
1922         IsUnionOfLinesPossible = Standard_False;
1923     }
1924     else
1925       IsUnionOfLinesPossible = Standard_False;
1926     if( c3d1->IsKind(STANDARD_TYPE(Geom_Circle)) && c3d2->IsKind(STANDARD_TYPE(Geom_Circle)) ) {
1927       Handle(Geom_Circle) C1 = Handle(Geom_Circle)::DownCast(c3d1);
1928       Handle(Geom_Circle) C2 = Handle(Geom_Circle)::DownCast(c3d2);
1929       gp_Pnt P01 = C1->Location();
1930       gp_Pnt P02 = C2->Location();
1931       if (P01.Distance(P02) > Precision::Confusion())
1932         IsUnionOfCirclesPossible = Standard_False;
1933     }
1934     else
1935       IsUnionOfCirclesPossible = Standard_False;
1936   }
1937   if (IsUnionOfLinesPossible && IsUnionOfCirclesPossible)
1938     return Standard_False;
1939
1940   //union of lines is possible
1941   if (IsUnionOfLinesPossible)
1942   {
1943     TopoDS_Vertex V[2];
1944     V[0] = sae.FirstVertex(TopoDS::Edge(theChain.First()));
1945     gp_Pnt PV1 = BRep_Tool::Pnt(V[0]);
1946     V[1] = sae.LastVertex(TopoDS::Edge(theChain.Last()));
1947     gp_Pnt PV2 = BRep_Tool::Pnt(V[1]);
1948     gp_Vec Vec(PV1, PV2);
1949     if (mySafeInputMode) {
1950       for (int k = 0; k < 2; k++) {
1951         if (!myContext->IsRecorded(V[k])) {
1952           TopoDS_Vertex Vcopy = TopoDS::Vertex(V[k].EmptyCopied());
1953           myContext->Replace(V[k], Vcopy);
1954           V[k] = Vcopy;
1955         }
1956         else
1957           V[k] = TopoDS::Vertex (myContext->Apply(V[k]));
1958       }
1959     }
1960     Handle(Geom_Line) L = new Geom_Line(gp_Ax1(PV1,Vec));
1961     Standard_Real dist = PV1.Distance(PV2);
1962     Handle(Geom_TrimmedCurve) tc = new Geom_TrimmedCurve(L,0.0,dist);
1963     TopoDS_Edge E;
1964     B.MakeEdge (E, tc ,Precision::Confusion());
1965     B.Add (E,V[0]);  B.Add (E,V[1]);
1966     B.UpdateVertex(V[0], 0., E, 0.);
1967     B.UpdateVertex(V[1], dist, E, 0.);
1968     UnionPCurves(theChain, E);
1969     OutEdge = E;
1970     return Standard_True;
1971   }
1972
1973   if (IsUnionOfCirclesPossible)
1974   {
1975     double f,l;
1976     TopoDS_Edge FE = TopoDS::Edge(theChain.First());
1977     Handle(Geom_Curve) c3d = BRep_Tool::Curve(FE,f,l);
1978
1979     if (c3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
1980       Handle(Geom_TrimmedCurve) tc =
1981         Handle(Geom_TrimmedCurve)::DownCast(c3d);
1982       c3d = tc->BasisCurve();
1983     }
1984     Handle(Geom_Circle) Cir = Handle(Geom_Circle)::DownCast(c3d);
1985
1986     TopoDS_Vertex V[2];
1987     V[0] = sae.FirstVertex(FE);
1988     V[1] = sae.LastVertex(TopoDS::Edge(theChain.Last()));
1989     TopoDS_Edge E;
1990     if (V[0].IsSame(V[1])) {
1991       // closed chain
1992       BRepAdaptor_Curve adef(FE);
1993       Handle(Geom_Circle) Cir1;
1994       double FP, LP;
1995       if ( FE.Orientation() == TopAbs_FORWARD)
1996       {
1997         FP = adef.FirstParameter();
1998         LP = adef.LastParameter();
1999       }
2000       else
2001       {
2002         FP = adef.LastParameter();
2003         LP = adef.FirstParameter();
2004       }
2005       if (Abs(FP) < Precision::PConfusion())
2006       {
2007         B.MakeEdge (E,Cir, Precision::Confusion());
2008         B.Add(E,V[0]);
2009         B.Add(E,V[1]);
2010         E.Orientation(FE.Orientation());
2011       }
2012       else
2013       {
2014         GC_MakeCircle MC1 (adef.Value(FP), adef.Value((FP + LP) * 0.5), adef.Value(LP));
2015         if (MC1.IsDone())
2016           Cir1 = MC1.Value();
2017         else
2018           return Standard_False;
2019         B.MakeEdge (E, Cir1, Precision::Confusion());
2020         B.Add(E,V[0]);
2021         B.Add(E,V[1]);
2022       }
2023     }
2024     else //open chain
2025     {
2026       Standard_Real ParamFirst = BRep_Tool::Parameter(V[0], FE);
2027       TopoDS_Vertex VertexLastOnFE = sae.LastVertex(FE);
2028       Standard_Real ParamLast  = BRep_Tool::Parameter(VertexLastOnFE, FE);
2029       
2030       if (mySafeInputMode) {
2031         for (int k = 0; k < 2; k++) {
2032           if (!myContext->IsRecorded(V[k])) {
2033             TopoDS_Vertex Vcopy = TopoDS::Vertex(V[k].EmptyCopied());
2034             myContext->Replace(V[k], Vcopy);
2035             V[k] = Vcopy;
2036           }
2037           else
2038             V[k] = TopoDS::Vertex (myContext->Apply(V[k]));
2039         }
2040       }
2041       
2042       gp_Pnt PointFirst = BRep_Tool::Pnt(V[0]);
2043       while (Abs(ParamLast - ParamFirst) > 7*M_PI/8)
2044         ParamLast = (ParamFirst + ParamLast)/2;
2045       BRepAdaptor_Curve BAcurveFE(FE);
2046       gp_Pnt PointLast = BAcurveFE.Value(ParamLast);
2047       gp_Pnt Origin = Cir->Circ().Location();
2048       gp_Dir Dir1 = gp_Vec(Origin, PointFirst);
2049       gp_Dir Dir2 = gp_Vec(Origin, PointLast);
2050       gp_Dir Vdir = Dir1 ^ Dir2;
2051       gp_Ax2 anAx2(Origin, Vdir, Dir1);
2052       Handle(Geom_Circle) aNewCircle = new Geom_Circle(anAx2, Cir->Radius());
2053       gp_Pnt PointLastInChain = BRep_Tool::Pnt(V[1]);
2054       gp_Dir DirLastInChain = gp_Vec(Origin, PointLastInChain);
2055       Standard_Real lpar = Dir1.AngleWithRef(DirLastInChain, Vdir);
2056       if (lpar < 0.)
2057         lpar += 2*M_PI;
2058
2059       Handle(Geom_TrimmedCurve) tc = new Geom_TrimmedCurve(aNewCircle,0.,lpar);
2060       B.MakeEdge (E,tc,Precision::Confusion());
2061       B.Add(E,V[0]);
2062       B.Add(E,V[1]);
2063       B.UpdateVertex(V[0], 0., E, 0.);
2064       B.UpdateVertex(V[1], lpar, E, 0.);
2065     }
2066     UnionPCurves(theChain, E);
2067     OutEdge = E;
2068     return Standard_True;
2069   }
2070   if (theChain.Length() > 1 && myConcatBSplines) {
2071     // second step: union edges with various curves
2072     // skl for bug 0020052 from Mantis: perform such unions
2073     // only if curves are bspline or bezier
2074
2075     TopoDS_Vertex VF = sae.FirstVertex(TopoDS::Edge(theChain.First()));
2076     TopoDS_Vertex VL = sae.LastVertex(TopoDS::Edge(theChain.Last()));
2077     Standard_Boolean NeedUnion = Standard_True;
2078     for(j = 1; j <= theChain.Length(); j++) {
2079       TopoDS_Edge edge = TopoDS::Edge(theChain.Value(j));
2080       TopLoc_Location Loc;
2081       Handle(Geom_Curve) c3d = BRep_Tool::Curve(edge,Loc,fp1,lp1);
2082       if(c3d.IsNull()) continue;
2083       if (c3d->IsKind(STANDARD_TYPE(Geom_TrimmedCurve))) {
2084         Handle(Geom_TrimmedCurve) tc =
2085           Handle(Geom_TrimmedCurve)::DownCast(c3d);
2086         c3d = tc->BasisCurve();
2087       }
2088       if( ( c3d->IsKind(STANDARD_TYPE(Geom_BSplineCurve)) ||
2089             c3d->IsKind(STANDARD_TYPE(Geom_BezierCurve)) ) ) continue;
2090       NeedUnion = Standard_False;
2091       break;
2092     }
2093     if(NeedUnion) {
2094 #ifdef OCCT_DEBUG
2095       std::cout<<"can not make analytical union => make approximation"<<std::endl;
2096 #endif
2097       TopoDS_Edge E = GlueEdgesWithPCurves(theChain, VF, VL);
2098       OutEdge = E;
2099       return Standard_True;
2100     }
2101     else {
2102 #ifdef OCCT_DEBUG
2103       std::cout<<"can not make approximation for such types of curves"<<std::endl;
2104 #endif
2105       return Standard_False;
2106     }
2107   }
2108   return Standard_False;
2109 }
2110
2111 //=======================================================================
2112 //function : IsMergingPossible
2113 //purpose  : Checks if merging of two edges is possible
2114 //=======================================================================
2115
2116 static Standard_Boolean IsMergingPossible(const TopoDS_Edge& edge1, const TopoDS_Edge& edge2, 
2117                                           double theAngTol, double theLinTol, 
2118                                           const TopTools_MapOfShape& AvoidEdgeVrt, const bool theLineDirectionOk,
2119                                           const gp_Pnt& theFirstPoint, const gp_Vec& theDirectionVec,
2120                                           const TopTools_IndexedDataMapOfShapeListOfShape& theVFmap)
2121 {
2122   Standard_Boolean IsDegE1 = BRep_Tool::Degenerated(edge1);
2123   Standard_Boolean IsDegE2 = BRep_Tool::Degenerated(edge2);
2124   
2125   if (IsDegE1 && IsDegE2)
2126   {
2127     //Find connstion point in 2d
2128     TopoDS_Face CommonFace;
2129     Standard_Real MinSqDist;
2130     TopAbs_Orientation OrOfE1OnFace, OrOfE2OnFace;
2131     Standard_Integer IndOnE1, IndOnE2;
2132     gp_Pnt2d PointsOnEdge1 [2], PointsOnEdge2 [2];
2133     if (!FindClosestPoints(edge1, edge2, theVFmap, CommonFace,
2134                            MinSqDist, OrOfE1OnFace, OrOfE2OnFace,
2135                            IndOnE1, IndOnE2, PointsOnEdge1, PointsOnEdge2))
2136       return Standard_False;
2137     
2138     if (MinSqDist <= Precision::SquareConfusion())
2139       return Standard_True;
2140     
2141     return Standard_False;
2142   }
2143   else if (IsDegE1 || IsDegE2)
2144     return Standard_False;
2145   
2146   TopoDS_Vertex CV = TopExp::LastVertex(edge1, Standard_True);
2147   if (CV.IsNull() || AvoidEdgeVrt.Contains(CV))
2148     return Standard_False;
2149
2150   BRepAdaptor_Curve ade1(edge1);
2151   BRepAdaptor_Curve ade2(edge2);
2152
2153   GeomAbs_CurveType t1 = ade1.GetType();
2154   GeomAbs_CurveType t2 = ade2.GetType();
2155
2156   if( t1 == GeomAbs_Circle && t2 == GeomAbs_Circle)
2157   {
2158     if (ade1.Circle().Location().Distance(ade2.Circle().Location()) > Precision::Confusion())
2159       return Standard_False;
2160   }
2161
2162   if( ( (t1 != GeomAbs_BezierCurve && t1 != GeomAbs_BSplineCurve) ||
2163       (t2 != GeomAbs_BezierCurve && t2 != GeomAbs_BSplineCurve)) && t1 != t2)
2164     return Standard_False;
2165
2166   gp_Vec Diff1, Diff2;
2167   gp_Pnt P1, P2;
2168   if (edge1.Orientation() == TopAbs_FORWARD)
2169     ade1.D1(ade1.LastParameter(), P1, Diff1);
2170   else
2171   {
2172     ade1.D1(ade1.FirstParameter(), P1, Diff1);
2173     Diff1 = -Diff1;
2174   }
2175
2176   if (edge2.Orientation() == TopAbs_FORWARD)
2177     ade2.D1(ade2.FirstParameter(), P2, Diff2);
2178   else
2179   {
2180     ade2.D1(ade2.LastParameter(), P2, Diff2);
2181     Diff2 = -Diff2;
2182   }
2183
2184   if (Diff1.Angle(Diff2) > theAngTol)
2185     return Standard_False;
2186
2187   if (theLineDirectionOk && t2 == GeomAbs_Line)
2188   {
2189     // Check that the accumulated deflection does not exceed the linear tolerance
2190     Standard_Real aLast = (edge2.Orientation() == TopAbs_FORWARD) ?
2191       ade2.LastParameter() : ade2.FirstParameter();
2192     gp_Vec aCurV(theFirstPoint, ade2.Value(aLast));
2193     Standard_Real aDD = theDirectionVec.CrossSquareMagnitude(aCurV);
2194     if (aDD > theLinTol*theLinTol)
2195       return Standard_False;
2196
2197     // Check that the accumulated angle does not exceed the angular tolerance.
2198     // For symmetry, check the angle between vectors of:
2199     // - first edge and resulting curve, and
2200     // - the last edge and resulting curve.
2201     if (theDirectionVec.Angle(aCurV) > theAngTol || Diff2.Angle(aCurV) > theAngTol)
2202       return Standard_False;
2203   }
2204
2205   return Standard_True;
2206 }
2207
2208 //=======================================================================
2209 //function : GetLineEdgePoints
2210 //purpose  : 
2211 //=======================================================================
2212 static Standard_Boolean GetLineEdgePoints(const TopoDS_Edge& theInpEdge, gp_Pnt& theFirstPoint, gp_Vec& theDirectionVec)
2213 {
2214   double f, l;
2215   Handle(Geom_Curve) aCur = BRep_Tool::Curve(theInpEdge, f, l);
2216   if(aCur.IsNull()) 
2217     return Standard_False;
2218
2219   Handle(Geom_TrimmedCurve) aTC = Handle(Geom_TrimmedCurve)::DownCast(aCur);
2220   if (!aTC.IsNull())
2221     aCur = aTC->BasisCurve();
2222
2223   if (aCur->DynamicType() != STANDARD_TYPE(Geom_Line))
2224     return Standard_False;
2225
2226   if (theInpEdge.Orientation() == TopAbs_REVERSED) {
2227     Standard_Real tmp = f;
2228     f = l;
2229     l = tmp;
2230   }
2231   theFirstPoint = aCur->Value(f);
2232   gp_Pnt aLP = aCur->Value(l);
2233   theDirectionVec = aLP.XYZ().Subtracted(theFirstPoint.XYZ());
2234   theDirectionVec.Normalize();
2235   return Standard_True;
2236 }
2237
2238 struct ShapeUpgrade_UnifySameDomain::SubSequenceOfEdges
2239 {
2240   TopTools_SequenceOfShape SeqsEdges;
2241   TopoDS_Edge UnionEdges;
2242 };
2243
2244 //=======================================================================
2245 //function : generateSubSeq
2246 //purpose  :
2247 //=======================================================================
2248 void ShapeUpgrade_UnifySameDomain::generateSubSeq (const TopTools_SequenceOfShape& anInpEdgeSeq,
2249                                                    NCollection_Sequence<SubSequenceOfEdges>& SeqOfSubSeqOfEdges,
2250                                                    Standard_Boolean IsClosed, double theAngTol, double theLinTol,
2251                                                    const TopTools_MapOfShape& AvoidEdgeVrt,
2252                                                    const TopTools_IndexedDataMapOfShapeListOfShape& theVFmap)
2253 {
2254   Standard_Boolean isOk = Standard_False;
2255   TopoDS_Edge edge1, edge2;
2256
2257   SubSequenceOfEdges SubSeq;
2258   TopoDS_Edge RefEdge = TopoDS::Edge(anInpEdgeSeq(1));
2259   SubSeq.SeqsEdges.Append(RefEdge);
2260   SeqOfSubSeqOfEdges.Append(SubSeq);
2261
2262   gp_Pnt aFirstPoint;
2263   gp_Vec aDirectionVec;
2264   Standard_Boolean isLineDirectionOk = GetLineEdgePoints(RefEdge, aFirstPoint, aDirectionVec);  
2265   
2266   for (int i = 1; i < anInpEdgeSeq.Length(); i++)
2267   {
2268     edge1 = TopoDS::Edge(anInpEdgeSeq(i));
2269     edge2 = TopoDS::Edge(anInpEdgeSeq(i+1));
2270     isOk = IsMergingPossible(edge1, edge2, theAngTol, theLinTol,
2271                              AvoidEdgeVrt, isLineDirectionOk, aFirstPoint, aDirectionVec, theVFmap);
2272     if (!isOk)
2273     {
2274       SubSequenceOfEdges aSubSeq;
2275       aSubSeq.SeqsEdges.Append(edge2);
2276       SeqOfSubSeqOfEdges.Append(aSubSeq);
2277       isLineDirectionOk = GetLineEdgePoints(edge2, aFirstPoint, aDirectionVec);
2278     }
2279     else
2280       SeqOfSubSeqOfEdges.ChangeLast().SeqsEdges.Append(edge2);
2281   }
2282   /// check first and last chain segments
2283   if (IsClosed && SeqOfSubSeqOfEdges.Length() > 1)
2284   {
2285     edge1 = TopoDS::Edge(anInpEdgeSeq.Last());
2286     edge2 = TopoDS::Edge(anInpEdgeSeq.First());
2287     if (IsMergingPossible(edge1, edge2, theAngTol, theLinTol,
2288                           AvoidEdgeVrt, Standard_False, aFirstPoint, aDirectionVec, theVFmap))
2289     {
2290       SeqOfSubSeqOfEdges.ChangeLast().SeqsEdges.Append(SeqOfSubSeqOfEdges.ChangeFirst().SeqsEdges);
2291       SeqOfSubSeqOfEdges.Remove(1);
2292     }
2293   }
2294 }
2295
2296 //=======================================================================
2297 //function : MergeEdges
2298 //purpose  : auxiliary
2299 //=======================================================================
2300 Standard_Boolean ShapeUpgrade_UnifySameDomain::MergeEdges(TopTools_SequenceOfShape& SeqEdges,
2301                                                           const TopTools_IndexedDataMapOfShapeListOfShape& theVFmap,
2302                                                           NCollection_Sequence<SubSequenceOfEdges>& SeqOfSubSeqOfEdges,
2303                                                           const TopTools_MapOfShape& NonMergVrt)
2304 {
2305   TopTools_IndexedDataMapOfShapeListOfShape aMapVE;
2306   Standard_Integer j;
2307   TopTools_MapOfShape VerticesToAvoid;
2308   const Standard_Integer aNbE = SeqEdges.Length();
2309   for (j = 1; j <= aNbE; j++)
2310   {
2311     TopoDS_Edge anEdge = TopoDS::Edge(SeqEdges(j));
2312     // fill in the map V-E
2313     for (TopoDS_Iterator it(anEdge.Oriented(TopAbs_FORWARD)); it.More(); it.Next())
2314     {
2315       TopoDS_Shape aV = it.Value();
2316       if (aV.Orientation() == TopAbs_FORWARD || aV.Orientation() == TopAbs_REVERSED)
2317       {
2318         if (!aMapVE.Contains(aV))
2319           aMapVE.Add(aV, TopTools_ListOfShape());
2320         aMapVE.ChangeFromKey(aV).Append(anEdge);
2321       }
2322     }
2323   }
2324   VerticesToAvoid.Unite(NonMergVrt);
2325
2326   // do loop while there are unused edges
2327   TopTools_MapOfShape aUsedEdges;
2328
2329   for (Standard_Integer iE = 1; iE <= aNbE; ++iE)
2330   {
2331     TopoDS_Edge edge = TopoDS::Edge (SeqEdges (iE));
2332     if (!aUsedEdges.Add (edge))
2333       continue;
2334
2335     // make chain for unite
2336     TopTools_SequenceOfShape aChain;
2337     aChain.Append(edge);
2338     TopoDS_Vertex V[2];
2339     TopExp::Vertices(edge, V[0], V[1], Standard_True);
2340
2341     // connect more edges to the chain in both directions
2342     for (j = 0; j < 2; j++)
2343     {
2344       Standard_Boolean isAdded = Standard_True;
2345       while (isAdded)
2346       {
2347         isAdded = Standard_False;
2348         if (V[j].IsNull())
2349           break;
2350         const TopTools_ListOfShape& aLE = aMapVE.FindFromKey(V[j]);
2351         for (TopTools_ListIteratorOfListOfShape itL(aLE); itL.More(); itL.Next())
2352         {
2353           edge = TopoDS::Edge(itL.Value());
2354           if (!aUsedEdges.Contains(edge))
2355           {
2356             TopoDS_Vertex V2[2];
2357             TopExp::Vertices(edge, V2[0], V2[1], Standard_True);
2358             // the neighboring edge must have V[j] reversed and located on the opposite end
2359             if (V2[1 - j].IsEqual(V[j].Reversed()))
2360             {
2361               if (j == 0)
2362                 aChain.Prepend(edge);
2363               else
2364                 aChain.Append(edge);
2365               aUsedEdges.Add(edge);
2366               V[j] = V2[j];
2367               isAdded = Standard_True;
2368               break;
2369             }
2370           }
2371         }
2372       }
2373     }
2374
2375     if (aChain.Length() < 2)
2376       continue;
2377
2378     Standard_Boolean IsClosed = Standard_False;
2379     if (V[0].IsSame ( V[1] ))
2380       IsClosed = Standard_True;
2381
2382     // split chain by vertices at which merging is not possible
2383     NCollection_Sequence<SubSequenceOfEdges> aOneSeq;
2384     generateSubSeq(aChain, aOneSeq, IsClosed, myAngTol, myLinTol, VerticesToAvoid, theVFmap);
2385
2386     // put sub-chains in the result
2387     SeqOfSubSeqOfEdges.Append(aOneSeq);
2388   }
2389
2390   for (int i = 1; i <= SeqOfSubSeqOfEdges.Length(); i++)
2391   {
2392     TopoDS_Edge UE;
2393     if (SeqOfSubSeqOfEdges(i).SeqsEdges.Length() < 2)
2394       continue;
2395     if (MergeSubSeq(SeqOfSubSeqOfEdges(i).SeqsEdges, theVFmap, UE))
2396       SeqOfSubSeqOfEdges(i).UnionEdges = UE;
2397   }
2398   return Standard_True;
2399 }
2400
2401 //=======================================================================
2402 //function : MergeSeq
2403 //purpose  : Tries to unify the sequence of edges with the set of
2404 //           another edges which lies on the same geometry
2405 //=======================================================================
2406 Standard_Boolean ShapeUpgrade_UnifySameDomain::MergeSeq (TopTools_SequenceOfShape& SeqEdges,
2407                                                          const TopTools_IndexedDataMapOfShapeListOfShape& theVFmap,
2408                                                          const TopTools_MapOfShape& nonMergVert)
2409 {
2410   NCollection_Sequence<SubSequenceOfEdges> SeqOfSubsSeqOfEdges;
2411   if (MergeEdges(SeqEdges, theVFmap, SeqOfSubsSeqOfEdges, nonMergVert))
2412   {
2413     for (Standard_Integer i = 1; i <= SeqOfSubsSeqOfEdges.Length(); i++ )
2414     {
2415       if (SeqOfSubsSeqOfEdges(i).UnionEdges.IsNull())
2416         continue;
2417
2418       myContext->Merge(SeqOfSubsSeqOfEdges(i).SeqsEdges,
2419         SeqOfSubsSeqOfEdges(i).UnionEdges);
2420     }
2421     return Standard_True;
2422   }
2423   return Standard_False;
2424 }
2425
2426 //=======================================================================
2427 //function : CheckSharedVertices
2428 //purpose  : Checks the sequence of edges on the presence of shared vertex 
2429 //=======================================================================
2430
2431 static void CheckSharedVertices(const TopTools_SequenceOfShape& theSeqEdges, 
2432                                 const TopTools_IndexedDataMapOfShapeListOfShape& theMapEdgesVertex,
2433                                 const TopTools_MapOfShape& theMapKeepShape,
2434                                 TopTools_MapOfShape& theShareVertMap)
2435 {
2436   ShapeAnalysis_Edge sae;
2437   TopTools_SequenceOfShape SeqVertexes;
2438   TopTools_MapOfShape MapVertexes;
2439   for (Standard_Integer k = 1; k <= theSeqEdges.Length(); k++ )
2440   {
2441     TopoDS_Vertex aV1 = sae.FirstVertex(TopoDS::Edge(theSeqEdges(k)));
2442     TopoDS_Vertex aV2 = sae.LastVertex(TopoDS::Edge(theSeqEdges(k)));
2443     if (!MapVertexes.Add(aV1))
2444       SeqVertexes.Append(aV1);
2445     if (!MapVertexes.Add(aV2))
2446       SeqVertexes.Append(aV2);
2447   }
2448
2449   for (Standard_Integer k = 1; k <= SeqVertexes.Length()/* && !IsSharedVertexPresent*/; k++ )
2450   {
2451     const TopTools_ListOfShape& ListEdgesV1 = theMapEdgesVertex.FindFromKey(SeqVertexes(k));
2452     if (ListEdgesV1.Extent() > 2 || theMapKeepShape.Contains(SeqVertexes(k)))
2453       theShareVertMap.Add(SeqVertexes(k));
2454   }
2455   //return theShareVertMap.IsEmpty() ? false : true;
2456 }
2457
2458 //=======================================================================
2459 //function : ShapeUpgrade_UnifySameDomain
2460 //purpose  : Constructor
2461 //=======================================================================
2462
2463 ShapeUpgrade_UnifySameDomain::ShapeUpgrade_UnifySameDomain()
2464   : myLinTol(Precision::Confusion()),
2465     myAngTol(Precision::Angular()),
2466     myUnifyFaces(Standard_True),
2467     myUnifyEdges (Standard_True),
2468     myConcatBSplines (Standard_False),
2469     myAllowInternal (Standard_False),
2470     mySafeInputMode(Standard_True),
2471     myHistory(new BRepTools_History)
2472 {
2473   myContext = new ShapeBuild_ReShape;
2474 }
2475
2476 //=======================================================================
2477 //function : ShapeUpgrade_UnifySameDomain
2478 //purpose  : Constructor
2479 //=======================================================================
2480
2481 ShapeUpgrade_UnifySameDomain::ShapeUpgrade_UnifySameDomain(const TopoDS_Shape& aShape,
2482                                                            const Standard_Boolean UnifyEdges,
2483                                                            const Standard_Boolean UnifyFaces,
2484                                                            const Standard_Boolean ConcatBSplines)
2485   : myInitShape (aShape),
2486     myLinTol(Precision::Confusion()),
2487     myAngTol(Precision::Angular()),
2488     myUnifyFaces(UnifyFaces),
2489     myUnifyEdges (UnifyEdges),
2490     myConcatBSplines (ConcatBSplines),
2491     myAllowInternal (Standard_False),
2492     mySafeInputMode (Standard_True),
2493     myShape (aShape),
2494     myHistory(new BRepTools_History)
2495 {
2496   myContext = new ShapeBuild_ReShape;
2497 }
2498
2499 //=======================================================================
2500 //function : Initialize
2501 //purpose  : 
2502 //=======================================================================
2503
2504 void ShapeUpgrade_UnifySameDomain::Initialize(const TopoDS_Shape& aShape,
2505                                               const Standard_Boolean UnifyEdges,
2506                                               const Standard_Boolean UnifyFaces,
2507                                               const Standard_Boolean ConcatBSplines)
2508 {
2509   myInitShape = aShape;
2510   myShape = aShape;
2511   myUnifyEdges = UnifyEdges;
2512   myUnifyFaces = UnifyFaces;
2513   myConcatBSplines = ConcatBSplines;
2514
2515   myContext->Clear();
2516   myKeepShapes.Clear();
2517   myFacePlaneMap.Clear();
2518   myEFmap.Clear();
2519   myFaceNewFace.Clear();
2520   myHistory->Clear();
2521 }
2522
2523 //=======================================================================
2524 //function : AllowInternalEdges
2525 //purpose  : 
2526 //=======================================================================
2527
2528 void ShapeUpgrade_UnifySameDomain::AllowInternalEdges (const Standard_Boolean theValue)
2529 {
2530   myAllowInternal = theValue;
2531 }
2532
2533 //=======================================================================
2534 //function : SetSafeInputMode
2535 //purpose  : 
2536 //=======================================================================
2537
2538 void ShapeUpgrade_UnifySameDomain::SetSafeInputMode(Standard_Boolean theValue)
2539 {
2540   mySafeInputMode = theValue;
2541 }
2542
2543 //=======================================================================
2544 //function : KeepShape
2545 //purpose  : 
2546 //=======================================================================
2547
2548 void ShapeUpgrade_UnifySameDomain::KeepShape(const TopoDS_Shape& theShape)
2549 {
2550   if (theShape.ShapeType() == TopAbs_EDGE || theShape.ShapeType() == TopAbs_VERTEX)
2551     myKeepShapes.Add(theShape);
2552 }
2553
2554 //=======================================================================
2555 //function : KeepShapes
2556 //purpose  : 
2557 //=======================================================================
2558
2559 void ShapeUpgrade_UnifySameDomain::KeepShapes(const TopTools_MapOfShape& theShapes)
2560 {
2561   for (TopTools_MapIteratorOfMapOfShape it(theShapes); it.More(); it.Next()) {
2562     if (it.Value().ShapeType() == TopAbs_EDGE || it.Value().ShapeType() == TopAbs_VERTEX)
2563       myKeepShapes.Add(it.Value());
2564   }
2565 }
2566
2567 //=======================================================================
2568 //function : UnifyFaces
2569 //purpose  : 
2570 //=======================================================================
2571
2572 void ShapeUpgrade_UnifySameDomain::UnifyFaces()
2573 {
2574   // creating map of edge faces for the whole shape
2575   TopTools_IndexedDataMapOfShapeListOfShape aGMapEdgeFaces;
2576   TopExp::MapShapesAndAncestors(myShape, TopAbs_EDGE, TopAbs_FACE, aGMapEdgeFaces);
2577   
2578   // unify faces in each shell separately
2579   TopExp_Explorer exps;
2580   for (exps.Init(myShape, TopAbs_SHELL); exps.More(); exps.Next())
2581     IntUnifyFaces(exps.Current(), aGMapEdgeFaces);
2582
2583   // gather all faces out of shells in one compound and unify them at once
2584   BRep_Builder aBB;
2585   TopoDS_Compound aCmp;
2586   aBB.MakeCompound(aCmp);
2587   Standard_Integer nbf = 0;
2588   for (exps.Init(myShape, TopAbs_FACE, TopAbs_SHELL); exps.More(); exps.Next(), nbf++)
2589     aBB.Add(aCmp, exps.Current());
2590
2591   if (nbf > 0)
2592     IntUnifyFaces(aCmp, aGMapEdgeFaces);
2593   
2594   myShape = myContext->Apply(myShape);
2595 }
2596
2597 //=======================================================================
2598 //function : SetFixWireModes
2599 //purpose  : 
2600 //=======================================================================
2601
2602 static void SetFixWireModes(ShapeFix_Face& theSff)
2603 {
2604   Handle(ShapeFix_Wire) aFixWire = theSff.FixWireTool();
2605   aFixWire->FixSelfIntersectionMode() = 0;
2606   aFixWire->FixNonAdjacentIntersectingEdgesMode() = 0;
2607   aFixWire->FixLackingMode() = 0;
2608   aFixWire->FixNotchedEdgesMode() = 0;
2609   aFixWire->ModifyTopologyMode() = Standard_False;
2610   aFixWire->ModifyRemoveLoopMode() = 0;
2611   aFixWire->FixGapsByRangesMode() = Standard_False;
2612   aFixWire->FixSmallMode() = 0;
2613 }
2614
2615 //=======================================================================
2616 //function : IntUnifyFaces
2617 //purpose  : 
2618 //=======================================================================
2619
2620 void ShapeUpgrade_UnifySameDomain::IntUnifyFaces(const TopoDS_Shape& theInpShape,
2621    TopTools_IndexedDataMapOfShapeListOfShape& theGMapEdgeFaces)
2622 {
2623   // creating map of edge faces for the shape
2624   TopTools_IndexedDataMapOfShapeListOfShape aMapEdgeFaces;
2625   TopExp::MapShapesAndAncestors(theInpShape, TopAbs_EDGE, TopAbs_FACE, aMapEdgeFaces);
2626
2627   // map of processed shapes
2628   TopTools_MapOfShape aProcessed;
2629
2630   // processing each face
2631   TopExp_Explorer exp;
2632   for (exp.Init(theInpShape, TopAbs_FACE); exp.More(); exp.Next()) {
2633     
2634     TopoDS_Face aFace = TopoDS::Face(exp.Current());
2635
2636     if (aProcessed.Contains(aFace))
2637       continue;
2638
2639     // Boundary edges for the new face
2640     TopTools_SequenceOfShape edges;
2641     TopTools_SequenceOfShape RemovedEdges;
2642
2643     Standard_Integer dummy;
2644     AddOrdinaryEdges(edges, aFace, dummy, RemovedEdges);
2645
2646     // Faces to get unified with the current faces
2647     TopTools_SequenceOfShape faces;
2648
2649     // Add the current face for unification
2650     faces.Append(aFace);
2651
2652     // surface and location to construct result
2653     TopLoc_Location aBaseLocation;
2654     Handle(Geom_Surface) aBaseSurface = BRep_Tool::Surface(aFace);
2655     aBaseSurface = ClearRts(aBaseSurface);
2656     TopAbs_Orientation RefFaceOrientation = aFace.Orientation();
2657
2658     //Take original surface
2659     TopoDS_Face RefFace;
2660     BRep_Builder BB;
2661     BB.MakeFace(RefFace, aBaseSurface, aBaseLocation, 0.);
2662     RefFace.Orientation(RefFaceOrientation);
2663     TopTools_MapOfShape MapEdgesWithTemporaryPCurves; //map of edges not lying on RefFace
2664     //these edges may be updated by temporary pcurves
2665     
2666     Standard_Real Uperiod = (aBaseSurface->IsUPeriodic())? aBaseSurface->UPeriod() : 0.;
2667     Standard_Real Vperiod = (aBaseSurface->IsVPeriodic())? aBaseSurface->VPeriod() : 0.;
2668
2669     // find adjacent faces to union
2670     Standard_Integer i;
2671     for (i = 1; i <= edges.Length(); i++) {
2672       TopoDS_Edge edge = TopoDS::Edge(edges(i));
2673       if (BRep_Tool::Degenerated(edge))
2674         continue;
2675
2676       // get connectivity of the edge in the global shape
2677       const TopTools_ListOfShape& aGList = theGMapEdgeFaces.FindFromKey(edge);
2678       if (!myAllowInternal && (aGList.Extent() != 2 || myKeepShapes.Contains(edge))) {
2679         // non manifold case is not processed unless myAllowInternal
2680         continue;
2681       }
2682       //
2683       // Get the faces connected through the edge in the current shape
2684       const TopTools_ListOfShape& aList = aMapEdgeFaces.FindFromKey(edge);
2685       if (aList.Extent() < 2) {
2686         continue;
2687       }
2688
2689       // for a planar face create and store pcurve of edge on face
2690       // to speed up all operations
2691       if (!mySafeInputMode && aBaseSurface->IsKind(STANDARD_TYPE(Geom_Plane)))
2692         BRepLib::BuildPCurveForEdgeOnPlane(edge, aFace);
2693
2694       // get normal of the face to compare it with normals of other faces
2695       gp_Dir aDN1;
2696       //
2697       // take intermediate point on edge to compute the normal
2698       Standard_Real f, l;
2699       BRep_Tool::Range(edge, f, l);
2700       Standard_Real aTMid = (f + l) * .5;
2701       //
2702       Standard_Boolean bCheckNormals = GetNormalToSurface(aFace, edge, aTMid, aDN1);
2703       //
2704       // Process the faces
2705       TopTools_ListIteratorOfListOfShape anIter(aList);
2706       for (; anIter.More(); anIter.Next()) {
2707         
2708         TopoDS_Face aCheckedFace = TopoDS::Face(anIter.Value());
2709         if (aCheckedFace.IsSame(aFace))
2710           continue;
2711
2712         if (aProcessed.Contains(aCheckedFace))
2713           continue;
2714
2715         if (bCheckNormals) {
2716           // get normal of checked face using the same parameter on edge
2717           gp_Dir aDN2;
2718           if (GetNormalToSurface(aCheckedFace, edge, aTMid, aDN2)) {
2719             // and check if the adjacent faces are having approximately same normals
2720             Standard_Real anAngle = aDN1.Angle(aDN2);
2721             if (anAngle > myAngTol) {
2722               continue;
2723             }
2724           }
2725         }
2726         //
2727         if (IsSameDomain(aFace,aCheckedFace, myLinTol, myAngTol, myFacePlaneMap)) {
2728
2729           if (AddOrdinaryEdges(edges, aCheckedFace, dummy, RemovedEdges)) {
2730             // sequence edges is modified
2731             i = dummy;
2732           }
2733
2734           faces.Append(aCheckedFace);
2735           aProcessed.Add(aCheckedFace);
2736           break;
2737         }
2738       }
2739     }
2740
2741     if (faces.Length() > 1) {
2742       if (myFacePlaneMap.IsBound(faces(1)))
2743       {
2744         const Handle(Geom_Plane)& aPlane = myFacePlaneMap(faces(1));
2745         TopLoc_Location aLoc;
2746         BB.UpdateFace(RefFace, aPlane, aLoc, Precision::Confusion());
2747       }
2748       //Add correct pcurves for the reference surface to the edges of other faces
2749       TopoDS_Face F_RefFace = RefFace;
2750       F_RefFace.Orientation(TopAbs_FORWARD);
2751       AddPCurves(faces, F_RefFace, MapEdgesWithTemporaryPCurves);
2752       
2753       // fill in the connectivity map for selected faces
2754       TopTools_IndexedDataMapOfShapeListOfShape aMapEF;
2755       for (i = 1; i <= faces.Length(); i++) {
2756         TopExp::MapShapesAndAncestors(faces(i), TopAbs_EDGE, TopAbs_FACE, aMapEF);
2757       }
2758       // Collect keep edges and multi-connected edges, i.e. edges that are internal to
2759       // the set of selected faces and have connections to other faces.
2760       TopTools_ListOfShape aKeepEdges;
2761       for (i = 1; i <= aMapEF.Extent(); i++) {
2762         const TopTools_ListOfShape& aLF = aMapEF(i);
2763         if (aLF.Extent() == 2) {
2764           const TopoDS_Shape& aE = aMapEF.FindKey(i);
2765           const TopTools_ListOfShape& aGLF = theGMapEdgeFaces.FindFromKey(aE);
2766           if (aGLF.Extent() > 2 || myKeepShapes.Contains(aE)) {
2767             aKeepEdges.Append(aE);
2768           }
2769         }
2770       } 
2771       if (!aKeepEdges.IsEmpty()) {
2772         if  (!myAllowInternal) {
2773           // Remove from the selection the faces which have no other connect edges 
2774           // and contain multi-connected edges and/or keep edges.
2775           TopTools_MapOfShape anAvoidFaces;
2776           TopTools_ListIteratorOfListOfShape it(aKeepEdges);
2777           for (; it.More(); it.Next()) {
2778             const TopoDS_Shape& aE = it.Value();
2779             const TopTools_ListOfShape& aLF = aMapEF.FindFromKey(aE);
2780             anAvoidFaces.Add(aLF.First());
2781             anAvoidFaces.Add(aLF.Last());
2782           }
2783           for (i = 1; i <= faces.Length(); i++) {
2784             if (anAvoidFaces.Contains(faces(i))) {
2785               // update the boundaries of merged area, for that
2786               // remove from 'edges' the edges of this face and add to 'edges' 
2787               // the edges of this face that were not present in 'edges' before
2788               Standard_Boolean hasConnectAnotherFaces = Standard_False;
2789               TopExp_Explorer ex(faces(i), TopAbs_EDGE);
2790               for (; ex.More() && !hasConnectAnotherFaces; ex.Next()) {
2791                 TopoDS_Shape aE = ex.Current();
2792                 const TopTools_ListOfShape& aLF = aMapEF.FindFromKey(aE);
2793                 if (aLF.Extent() > 1) {
2794                   for (it.Init(aLF); it.More() && !hasConnectAnotherFaces; it.Next()) {
2795                     if (!anAvoidFaces.Contains(it.Value()))
2796                       hasConnectAnotherFaces = Standard_True;
2797                   }
2798                 }
2799               }
2800               if (!hasConnectAnotherFaces) {
2801                 AddOrdinaryEdges(edges, faces(i), dummy, RemovedEdges);
2802                 faces.Remove(i);
2803                 i--;
2804               }
2805             }
2806           }
2807           // check if the faces with keep edges contained in 
2808           // already updated the boundaries of merged area
2809           if (!faces.IsEmpty()) {
2810             TopTools_MapOfShape aMapFaces;
2811             for (i = 1; i <= faces.Length(); i++) {
2812               aMapFaces.Add(faces(i));
2813             }
2814             for (it.Init(aKeepEdges); it.More(); it.Next()) {
2815               const TopoDS_Shape& aE = it.Value();
2816               const TopTools_ListOfShape& aLF = aMapEF.FindFromKey(aE);
2817               if (aLF.Extent() < 2)
2818                 continue;
2819               if (aMapFaces.Contains(aLF.First()) && 
2820                   aMapFaces.Contains(aLF.Last())) {
2821                 for (i = 1; i <= faces.Length(); i++) {
2822                   if (faces(i).IsEqual(aLF.First()) ||
2823                       faces(i).IsEqual(aLF.Last())) {
2824                     AddOrdinaryEdges(edges, faces(i), dummy, RemovedEdges);
2825                     faces.Remove(i);
2826                     i--;
2827                   }
2828                 }
2829               }
2830             }
2831           }
2832         } //if (!myAllowInternal)
2833         else { //internal edges are allowed
2834           // add multi-connected and keep edges as internal in new face
2835           TopTools_ListIteratorOfListOfShape it(aKeepEdges);
2836           for (; it.More(); it.Next()) {
2837             const TopoDS_Shape& aE = it.Value();
2838             edges.Append(aE.Oriented(TopAbs_INTERNAL));
2839           }
2840         }
2841       } //if (!aKeepEdges.IsEmpty())
2842     } //if (faces.Length() > 1)
2843
2844     TopTools_IndexedDataMapOfShapeListOfShape aMapEF;
2845     for (i = 1; i <= faces.Length(); i++)
2846       TopExp::MapShapesAndUniqueAncestors(faces(i), TopAbs_EDGE, TopAbs_FACE, aMapEF);
2847     
2848     //Correct orientation of edges
2849     for (Standard_Integer ii = 1; ii <= edges.Length(); ii++)
2850     {
2851       const TopoDS_Shape& anEdge = edges (ii);
2852       Standard_Integer indE = aMapEF.FindIndex (anEdge);
2853       const TopTools_ListOfShape& aLF = aMapEF (indE);
2854       if (myAllowInternal &&
2855           myKeepShapes.Contains(anEdge) &&
2856           aLF.Extent() == 2)
2857         edges(ii).Orientation(TopAbs_INTERNAL);
2858       
2859       if (anEdge.Orientation() != TopAbs_INTERNAL)
2860         edges (ii) = aMapEF.FindKey (indE);
2861     }
2862
2863     //Exclude internal edges
2864     TopTools_IndexedMapOfShape InternalEdges;
2865     Standard_Integer ind_e = 1;
2866     while (ind_e <= edges.Length())
2867     {
2868       const TopoDS_Shape& anEdge = edges(ind_e);
2869       if (anEdge.Orientation() == TopAbs_INTERNAL)
2870       {
2871         InternalEdges.Add(anEdge);
2872         edges.Remove(ind_e);
2873       }
2874       else
2875         ind_e++;
2876     }    
2877
2878     if (RefFaceOrientation == TopAbs_REVERSED)
2879       for (Standard_Integer ii = 1; ii <= edges.Length(); ii++)
2880         edges(ii).Reverse();
2881     TopoDS_Face F_RefFace = RefFace;
2882     F_RefFace.Orientation(TopAbs_FORWARD);
2883     
2884     // all faces collected in the sequence. Perform union of faces
2885     if (faces.Length() > 1)
2886     {
2887       Standard_Real CoordTol = Precision::Confusion();
2888       TopTools_MapOfShape edgesMap;
2889       CoordTol = ComputeMinEdgeSize(edges, F_RefFace, edgesMap);
2890       CoordTol /= 10.;
2891       CoordTol = Max(CoordTol, Precision::Confusion());
2892
2893       TopTools_IndexedDataMapOfShapeListOfShape VEmap;
2894       for (Standard_Integer ind = 1; ind <= edges.Length(); ind++)
2895         TopExp::MapShapesAndUniqueAncestors(edges(ind), TopAbs_VERTEX, TopAbs_EDGE, VEmap);
2896       
2897       //Perform relocating to new U-origin
2898       //Define boundaries in 2d space of RefFace
2899       if (Uperiod != 0.)
2900       {
2901         //try to find a real seam edge - if it exists, do nothing
2902         Standard_Boolean SeamFound = Standard_False;
2903         for (Standard_Integer ii = 1; ii <= faces.Length(); ii++)
2904         {
2905           const TopoDS_Face& face_ii = TopoDS::Face(faces(ii));
2906           TopoDS_Wire anOuterWire = BRepTools::OuterWire(face_ii);
2907           TopoDS_Iterator itw(anOuterWire);
2908           for (; itw.More(); itw.Next())
2909           {
2910             const TopoDS_Edge& anEdge = TopoDS::Edge(itw.Value());
2911             if (BRepTools::IsReallyClosed(anEdge, face_ii))
2912             {
2913               SeamFound = Standard_True;
2914               break;
2915             }
2916           }
2917         }
2918         
2919         if (!SeamFound)
2920         {
2921           //try to find the origin of U in 2d space
2922           //so that all the faces are in [origin, origin + Uperiod]
2923           Standard_Real Umax;
2924           Standard_Integer i_face_max;
2925           FindFaceWithMaxUbound(faces, F_RefFace, edgesMap, Umax, i_face_max);
2926           
2927           TopTools_MapOfShape UsedEdges;
2928           NCollection_DataMap<TopoDS_Shape, Handle(Geom2d_Curve)> EdgeNewPCurve;
2929
2930           //Relocate pcurves to new U-origin
2931           RelocatePCurvesToNewUorigin(edges, faces(i_face_max), F_RefFace, CoordTol, Uperiod,
2932                                       VEmap, EdgeNewPCurve, UsedEdges);
2933
2934           //PCurves from unused edges (may be degenerated edges)
2935           for (Standard_Integer ind = 1; ind <= edges.Length(); ind++)
2936           {
2937             const TopoDS_Edge& anEdge = TopoDS::Edge(edges(ind));
2938             if (!UsedEdges.Contains(anEdge))
2939             {
2940               Standard_Real fpar, lpar;
2941               Handle(Geom2d_Curve) aPCurve = BRep_Tool::CurveOnSurface(anEdge, F_RefFace, fpar, lpar);
2942               aPCurve = new Geom2d_TrimmedCurve(aPCurve, fpar, lpar);
2943               EdgeNewPCurve.Bind(anEdge, aPCurve);
2944             }
2945           }
2946           
2947           //Restore VEmap
2948           VEmap.Clear();
2949           for (Standard_Integer ind = 1; ind <= edges.Length(); ind++)
2950             TopExp::MapShapesAndUniqueAncestors(edges(ind), TopAbs_VERTEX, TopAbs_EDGE, VEmap);
2951
2952           //Find NewUmin and NewUmax
2953           Standard_Real NewUmin = RealLast(), NewUmax = RealFirst();
2954           for (Standard_Integer ii = 1; ii <= edges.Length(); ii++)
2955           {
2956             const Handle(Geom2d_Curve)& aPCurve = EdgeNewPCurve(edges(ii));
2957             UpdateBoundaries(aPCurve, aPCurve->FirstParameter(), aPCurve->LastParameter(),
2958                              NewUmin, NewUmax);
2959           }
2960           
2961           if (NewUmax - NewUmin < Uperiod - CoordTol &&
2962               !(-Precision::Confusion() < NewUmin && NewUmin < Uperiod+Precision::Confusion() &&
2963                 -Precision::Confusion() < NewUmax && NewUmax < Uperiod+Precision::Confusion()))
2964           {
2965             //we can build a face without seam edge:
2966             //update the edges with earlier computed relocated pcurves
2967             //fitting into (NewUorigin, NewUorigin + Uperiod)
2968             Standard_Real umin, umax, vmin, vmax;
2969             aBaseSurface->Bounds(umin, umax, vmin, vmax);
2970             Standard_Real RestSpaceInU = Uperiod - (NewUmax - NewUmin);
2971             Standard_Real NewUorigin = NewUmin - RestSpaceInU/2;
2972             if (NewUorigin < umin)
2973               NewUorigin = umin;
2974             Handle(Geom_Surface) NewSurf;
2975             if (NewUorigin == umin)
2976               NewSurf = aBaseSurface;
2977             else
2978               NewSurf = new Geom_RectangularTrimmedSurface(aBaseSurface,
2979                                                            NewUorigin, NewUorigin + Uperiod,
2980                                                            Standard_True); //trim U
2981             TopoDS_Face OldRefFace = RefFace;
2982             Handle(Geom2d_Curve) NullPCurve;
2983             RefFace.Nullify();
2984             BB.MakeFace(RefFace, NewSurf, aBaseLocation, 0.);
2985             for (Standard_Integer ii = 1; ii <= edges.Length(); ii++)
2986             {
2987               TopoDS_Edge anEdge = TopoDS::Edge(edges(ii));
2988               if (MapEdgesWithTemporaryPCurves.Contains(anEdge))
2989                 BB.UpdateEdge(anEdge, NullPCurve, OldRefFace, 0.);
2990               const Handle(Geom2d_Curve)& aPCurve = EdgeNewPCurve(anEdge);
2991               BB.UpdateEdge(anEdge, aPCurve, RefFace, 0.);
2992             }
2993           }
2994         } //if (!SeamFound)
2995       } //if (Uperiod != 0.)
2996       ////////////////////////////////////
2997       F_RefFace = RefFace;
2998       F_RefFace.Orientation(TopAbs_FORWARD);
2999       
3000       TopTools_SequenceOfShape NewFaces, NewWires;
3001       
3002       if (Uperiod == 0 || Vperiod == 0)
3003       {
3004         //Set the "periods" for closed non-periodic surface
3005         TopLoc_Location aLoc;
3006         Handle(Geom_Surface) aSurf = BRep_Tool::Surface(RefFace, aLoc);
3007         if (aSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface)))
3008           aSurf = (Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurf))->BasisSurface();
3009         Standard_Real Ufirst, Ulast, Vfirst, Vlast;
3010         aSurf->Bounds(Ufirst, Ulast, Vfirst, Vlast);
3011         if (Uperiod == 0 && aSurf->IsUClosed())
3012           Uperiod = Ulast - Ufirst;
3013         if (Vperiod == 0 && aSurf->IsVClosed())
3014           Vperiod = Vlast - Vfirst;
3015       }
3016
3017       TopTools_MapOfShape UsedEdges;
3018
3019       Standard_Real FaceUmin = RealLast();
3020       Standard_Real FaceVmin = RealLast();
3021       for (Standard_Integer ii = 1; ii <= edges.Length(); ii++)
3022       {
3023         const TopoDS_Edge& anEdge = TopoDS::Edge(edges(ii));
3024         BRepAdaptor_Curve2d aBAcurve(anEdge, F_RefFace);
3025         gp_Pnt2d aFirstPoint = aBAcurve.Value(aBAcurve.FirstParameter());
3026         gp_Pnt2d aLastPoint  = aBAcurve.Value(aBAcurve.LastParameter());
3027         
3028         if (aFirstPoint.X() < FaceUmin)
3029           FaceUmin = aFirstPoint.X();
3030         if (aLastPoint.X() < FaceUmin)
3031           FaceUmin = aLastPoint.X();
3032         
3033         if (aFirstPoint.Y() < FaceVmin)
3034           FaceVmin = aFirstPoint.Y();
3035         if (aLastPoint.Y() < FaceVmin)
3036           FaceVmin = aLastPoint.Y();
3037       }
3038
3039       //Building new wires from <edges>
3040       //and build faces
3041       while (!edges.IsEmpty())
3042       {
3043         //try to find non-degenerated edge
3044         TopoDS_Edge StartEdge = TopoDS::Edge(edges(1));
3045         Standard_Integer istart = 1;
3046         while (BRep_Tool::Degenerated(StartEdge) &&
3047                istart < edges.Length())
3048         {
3049           istart++;
3050           StartEdge = TopoDS::Edge(edges(istart));
3051         }
3052
3053         TopoDS_Wire aNewWire;
3054         BB.MakeWire(aNewWire);
3055         BB.Add(aNewWire, StartEdge);
3056         RemoveEdgeFromMap(StartEdge, VEmap);
3057         TopTools_IndexedMapOfShape SplittingVertices;
3058         
3059         Standard_Real fpar, lpar;
3060         Handle(Geom2d_Curve) StartPCurve = BRep_Tool::CurveOnSurface(StartEdge, F_RefFace, fpar, lpar);
3061         TopoDS_Vertex StartVertex, CurVertex;
3062         TopExp::Vertices(StartEdge, StartVertex, CurVertex, Standard_True); //with orientation
3063         Standard_Real StartParam, CurParam;
3064         if (StartEdge.Orientation() == TopAbs_FORWARD)
3065         {
3066           StartParam = fpar; CurParam = lpar;
3067         }
3068         else
3069         {
3070           StartParam = lpar; CurParam = fpar;
3071         }
3072         gp_Pnt2d StartPoint = StartPCurve->Value(StartParam);
3073         gp_Pnt2d CurPoint   = StartPCurve->Value(CurParam);
3074         
3075         TopoDS_Edge CurEdge = StartEdge;
3076         for (;;) //loop till the end of current new wire
3077         {
3078           TopoDS_Edge NextEdge;
3079           gp_Pnt2d NextPoint;
3080           
3081           const TopTools_ListOfShape& Elist = VEmap.FindFromKey(CurVertex);
3082           TopTools_ListIteratorOfListOfShape itl(Elist);
3083           if (Elist.IsEmpty())
3084           {
3085             if (CurVertex.IsSame(StartVertex))
3086             {
3087               //Points of two vertices coincide in 3d but may be not in 2d
3088               if (Uperiod != 0. &&
3089                   Abs(StartPoint.X() - CurPoint.X()) > Uperiod/2) //end of parametric space
3090               {
3091                 //<edges> do not contain seams => we must reconstruct the seam up to <NextEdge>
3092                 gp_Pnt2d StartOfNextEdge;
3093                 TopoDS_Vertex LastVertexOfSeam;
3094                 ReconstructMissedSeam(edges, RemovedEdges, UsedEdges, F_RefFace, CurVertex,
3095                                       CurPoint, Uperiod, FaceUmin, Standard_True, CoordTol,
3096                                       NextEdge, aNewWire, NextPoint,
3097                                       StartOfNextEdge, LastVertexOfSeam, VEmap);
3098               }
3099               else if (Vperiod != 0. &&
3100                        Abs(StartPoint.Y() - CurPoint.Y()) > Vperiod/2) //end of parametric space
3101               {
3102                 //<edges> do not contain seams => we must reconstruct the seam up to <NextEdge>
3103                 gp_Pnt2d StartOfNextEdge;
3104                 TopoDS_Vertex LastVertexOfSeam;
3105                 ReconstructMissedSeam(edges, RemovedEdges, UsedEdges, F_RefFace, CurVertex,
3106                                       CurPoint, Vperiod, FaceVmin, Standard_False, CoordTol,
3107                                       NextEdge, aNewWire, NextPoint,
3108                                       StartOfNextEdge, LastVertexOfSeam, VEmap);
3109               }
3110               else
3111               {
3112                 break; //end of wire
3113               }
3114             }
3115           }
3116           
3117           if (NextEdge.IsNull())
3118           {
3119             Standard_Boolean EndOfWire = Standard_False;
3120             
3121             TopTools_ListOfShape TmpElist, TrueElist;
3122             //<TrueElist> will be the list of candidates to become <NextEdge>
3123             for (itl.Initialize(Elist); itl.More(); itl.Next())
3124             {
3125               const TopoDS_Edge& anEdge = TopoDS::Edge(itl.Value());
3126               if (UsedEdges.Contains(anEdge))
3127                 continue;
3128               TopoDS_Vertex aFirstVertex = TopExp::FirstVertex(anEdge, Standard_True);
3129               if (!aFirstVertex.IsSame(CurVertex))
3130                 continue;
3131               TmpElist.Append(anEdge);
3132             }
3133             if (TmpElist.Extent() <= 1 ||
3134                 (Uperiod != 0. || Vperiod != 0))
3135               TrueElist.Assign(TmpElist);
3136             else
3137             {
3138               //we must choose the closest direction - the biggest angle
3139               SplittingVertices.Add (CurVertex);
3140               Standard_Real MaxAngle = RealFirst();
3141               TopoDS_Edge TrueEdge;
3142               Handle(Geom2d_Curve) CurPCurve = BRep_Tool::CurveOnSurface(CurEdge, F_RefFace, fpar, lpar);
3143               CurParam = (CurEdge.Orientation() == TopAbs_FORWARD)? lpar : fpar;
3144               gp_Vec2d CurDir;
3145               CurPCurve->D1(CurParam, CurPoint, CurDir);
3146               CurDir.Normalize();
3147               if (CurEdge.Orientation() == TopAbs_REVERSED)
3148                 CurDir.Reverse();
3149               for (itl.Initialize(TmpElist); itl.More(); itl.Next())
3150               {
3151                 const TopoDS_Edge& anEdge = TopoDS::Edge(itl.Value());
3152                 Handle(Geom2d_Curve) aPCurve = BRep_Tool::CurveOnSurface(anEdge, F_RefFace, fpar, lpar);
3153                 Standard_Real aParam = (anEdge.Orientation() == TopAbs_FORWARD)? fpar : lpar;
3154                 gp_Pnt2d aPoint;
3155                 gp_Vec2d aDir;
3156                 aPCurve->D1(aParam, aPoint, aDir);
3157                 aDir.Normalize();
3158                 if (anEdge.Orientation() == TopAbs_REVERSED)
3159                   aDir.Reverse();
3160                 Standard_Real anAngle = CurDir.Angle(aDir);
3161                 if (anAngle > MaxAngle)
3162                 {
3163                   MaxAngle = anAngle;
3164                   TrueEdge = anEdge;
3165                 }
3166               }
3167               TrueElist.Append(TrueEdge);
3168             }
3169
3170             //Find next edge in TrueElist
3171             for (itl.Initialize(TrueElist); itl.More(); itl.Next())
3172             {
3173               const TopoDS_Edge& anEdge = TopoDS::Edge(itl.Value());
3174               
3175               Handle(Geom2d_Curve) aPCurve = BRep_Tool::CurveOnSurface(anEdge, F_RefFace, fpar, lpar);
3176               Standard_Real aParam = (anEdge.Orientation() == TopAbs_FORWARD)? fpar : lpar;
3177               gp_Pnt2d aPoint = aPCurve->Value(aParam);
3178               Standard_Real DiffU = Abs(aPoint.X() - CurPoint.X());
3179               Standard_Real DiffV = Abs(aPoint.Y() - CurPoint.Y());
3180               if (Uperiod != 0. &&
3181                   DiffU > CoordTol &&
3182                   Abs(DiffU - Uperiod) > CoordTol) //may be it is a deg.vertex
3183                 continue;
3184               if (Vperiod != 0. &&
3185                   DiffV > CoordTol &&
3186                   Abs(DiffV - Vperiod) > CoordTol) //may be it is a deg.vertex
3187                 continue;
3188               
3189               //Check: may be <CurPoint> and <aPoint> are on Period from each other
3190               if (Uperiod != 0. && DiffU > Uperiod/2) //end of parametric space
3191               {
3192                 //<edges> do not contain seams => we must reconstruct the seam up to <NextEdge>
3193                 gp_Pnt2d StartOfNextEdge;
3194                 TopoDS_Vertex LastVertexOfSeam;
3195                 ReconstructMissedSeam(edges, RemovedEdges, UsedEdges, F_RefFace, CurVertex,
3196                                       CurPoint, Uperiod, FaceUmin, Standard_True, CoordTol,
3197                                       NextEdge, aNewWire, NextPoint,
3198                                       StartOfNextEdge, LastVertexOfSeam, VEmap);
3199                 
3200                 //Check: may be it is the end
3201                 if (LastVertexOfSeam.IsSame(StartVertex) &&
3202                     Abs(StartPoint.X() - StartOfNextEdge.X()) < Uperiod/2)
3203                   EndOfWire = Standard_True;
3204                 
3205                 break;
3206               } //if (Uperiod != 0. && Abs(aPoint.X() - CurPoint.X()) > Uperiod/2)
3207               else if (Vperiod != 0. && DiffV > Vperiod/2) //end of parametric space
3208               {
3209                 //<edges> do not contain seams => we must reconstruct the seam up to <NextEdge>
3210                 gp_Pnt2d StartOfNextEdge;
3211                 TopoDS_Vertex LastVertexOfSeam;
3212                 ReconstructMissedSeam(edges, RemovedEdges, UsedEdges, F_RefFace, CurVertex,
3213                                       CurPoint, Vperiod, FaceVmin, Standard_False, CoordTol,
3214                                       NextEdge, aNewWire, NextPoint,
3215                                       StartOfNextEdge, LastVertexOfSeam, VEmap);
3216                 
3217                 //Check: may be it is the end
3218                 if (LastVertexOfSeam.IsSame(StartVertex) &&
3219                     Abs(StartPoint.Y() - StartOfNextEdge.Y()) < Vperiod/2)
3220                   EndOfWire = Standard_True;
3221                 
3222                 break;
3223               }
3224               else
3225               {
3226                 NextEdge = anEdge;
3227                 Standard_Real LastParam = (NextEdge.Orientation() == TopAbs_FORWARD)? lpar : fpar;
3228                 NextPoint = aPCurve->Value(LastParam);
3229                 break;
3230               }
3231             } //for (itl.Initialize(TrueElist); itl.More(); itl.Next())
3232             
3233             if (EndOfWire)
3234               break;
3235           }
3236           
3237           if (NextEdge.IsNull())
3238           {
3239             //throw Standard_ConstructionError("Construction of unified wire failed: no proper edge found");
3240             return;
3241           }
3242           else
3243           {
3244             CurPoint = NextPoint;
3245             CurEdge = NextEdge;
3246             CurVertex = TopExp::LastVertex(CurEdge, Standard_True); //with orientation
3247             BB.Add(aNewWire, CurEdge);
3248             UsedEdges.Add(CurEdge);
3249             RemoveEdgeFromMap(CurEdge, VEmap);
3250           }
3251         } //for (;;)
3252         aNewWire.Closed(Standard_True);
3253         UsedEdges.Add(StartEdge);
3254         
3255         //Remove used edges from sequence
3256         Standard_Integer ind = 1;