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