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