0024059: Eliminate compiler warning C4701 in MSVC++ with warning level 4
[occt.git] / src / BOPTools / BOPTools_AlgoTools.cxx
1 // Created by: Peter KURNEV
2 // Copyright (c) 2010-2012 OPEN CASCADE SAS
3 // Copyright (c) 2007-2010 CEA/DEN, EDF R&D, OPEN CASCADE
4 // Copyright (c) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, CEDRAT,
5 //                         EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
6 //
7 // The content of this file is subject to the Open CASCADE Technology Public
8 // License Version 6.5 (the "License"). You may not use the content of this file
9 // except in compliance with the License. Please obtain a copy of the License
10 // at http://www.opencascade.org and read it completely before using this file.
11 //
12 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
13 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
14 //
15 // The Original Code and all software distributed under the License is
16 // distributed on an "AS IS" basis, without warranty of any kind, and the
17 // Initial Developer hereby disclaims all such warranties, including without
18 // limitation, any warranties of merchantability, fitness for a particular
19 // purpose or non-infringement. Please see the License for the specific terms
20 // and conditions governing the rights and limitations under the License.
21
22
23 #include <BOPTools_AlgoTools.ixx>
24 //
25 #include <Precision.hxx>
26 //
27 #include <gp_Pnt.hxx>
28 #include <gp_XYZ.hxx>
29 #include <gp_Pnt2d.hxx>
30 #include <gp_Cylinder.hxx>
31 #include <gp_Cone.hxx>
32 #include <gp_Sphere.hxx>
33 #include <gp_Torus.hxx>
34 //
35 #include <Geom2d_Curve.hxx>
36 #include <Geom_Surface.hxx>
37 #include <Geom_Plane.hxx>
38 #include <Geom_TrimmedCurve.hxx>
39 #include <Geom_Curve.hxx>
40 #include <GeomAPI_ProjectPointOnSurf.hxx>
41 #include <Geom2dInt_Geom2dCurveTool.hxx>
42 #include <GeomAdaptor_Surface.hxx>
43 //
44 #include <TopAbs_Orientation.hxx>
45 //
46 #include <TopoDS_Compound.hxx>
47 #include <TopoDS_CompSolid.hxx>
48 #include <TopoDS_Solid.hxx>
49 #include <TopoDS_Shell.hxx>
50 #include <TopoDS_Wire.hxx>
51 //
52 #include <BRep_Builder.hxx>
53 #include <BRep_Tool.hxx>
54 #include <BRepLib.hxx>
55 #include <BRepAdaptor_Curve2d.hxx>
56 #include <BRepClass3d_SolidClassifier.hxx>
57 #include <TopExp.hxx>
58 #include <TopExp_Explorer.hxx>
59 //
60 #include <IntTools_Tools.hxx>
61 //
62 #include <BOPTools.hxx>
63 #include <BOPTools_CoupleOfShape.hxx>
64 #include <BOPTools_ListOfCoupleOfShape.hxx>
65 #include <BOPTools_AlgoTools2D.hxx>
66 #include <BOPTools_AlgoTools3D.hxx>
67 //
68 #include <BOPCol_IndexedMapOfShape.hxx>
69 #include <BOPCol_MapOfShape.hxx>
70 //
71 #include <BOPInt_ShrunkRange.hxx>
72 //
73
74 static
75   Standard_Real AngleWithRef(const gp_Dir& theD1,
76                              const gp_Dir& theD2,
77                              const gp_Dir& theDRef);
78 static 
79   inline Standard_Real fsqrt(Standard_Real val);
80
81 static
82   Standard_Boolean FindFacePairs (const TopoDS_Edge& theE,
83                                   const BOPCol_ListOfShape& thLF,
84                                   BOPTools_ListOfCoupleOfShape& theLCFF,
85                                   Handle(BOPInt_Context)& theContext);
86 static
87   TopAbs_Orientation Orientation(const TopoDS_Edge& anE,
88                                  const TopoDS_Face& aF);
89
90 static 
91   void GetFaceDir(const TopoDS_Edge& aE,
92                   const TopoDS_Face& aF,
93                   const gp_Pnt& aP,
94                   const Standard_Real aT,
95                   const gp_Dir& aDTgt,
96                   gp_Dir& aDN,
97                   gp_Dir& aDB,
98                   Handle(BOPInt_Context)& theContext);
99 static 
100   Standard_Boolean FindPointInFace(const TopoDS_Edge& aE,
101                                    const TopoDS_Face& aF,
102                                    const gp_Pnt& aP,
103                                    const Standard_Real aT,
104                                    gp_Dir& aDB,
105                                    gp_Pnt& aPOut,
106                                    Handle(BOPInt_Context)& theContext,
107                                    const GeomAdaptor_Surface& aGAS);
108
109 //=======================================================================
110 // function: MakeConnexityBlocks
111 // purpose: 
112 //=======================================================================
113   void BOPTools_AlgoTools::MakeConnexityBlocks (const TopoDS_Shape& theS,
114                                             const TopAbs_ShapeEnum theType1,
115                                             const TopAbs_ShapeEnum theType2,
116                                             BOPCol_ListOfShape& theLCB)
117 {
118   Standard_Integer  aNbF, aNbAdd, aNbAdd1, i;
119   BRep_Builder aBB;
120   TopoDS_Compound aC;
121   TopoDS_Iterator aIt;
122   TopExp_Explorer aExp;
123   BOPCol_MapOfShape aMP;
124   BOPCol_IndexedMapOfShape aMCB, aMAdd, aMAdd1;
125   BOPCol_IndexedDataMapOfShapeListOfShape aMEF;
126   BOPCol_ListIteratorOfListOfShape aItLF;
127   //
128   // 1. aMEF
129   BOPTools::MapShapesAndAncestors(theS, theType1, theType2, aMEF);
130   //
131   // 2. aMCB
132   aIt.Initialize(theS);
133   for (; aIt.More(); aIt.Next()) {
134     const TopoDS_Shape& aF1=aIt.Value();
135     if (aMP.Contains(aF1)) {
136       continue;
137     }
138     //
139     aMCB.Clear();
140     aMAdd.Clear();
141     aMAdd.Add(aF1);
142     //
143     for(;;) {
144       aMAdd1.Clear();
145       //
146       aNbAdd = aMAdd.Extent();
147       for (i=1; i<=aNbAdd; ++i) {
148         const TopoDS_Shape& aF=aMAdd(i);
149         //
150         aExp.Init(aF, theType1);
151         for (; aExp.More(); aExp.Next()) {
152           const TopoDS_Shape& aE=aExp.Current();
153           //
154           const BOPCol_ListOfShape& aLF=aMEF.FindFromKey(aE);
155           aItLF.Initialize(aLF);
156           for (; aItLF.More(); aItLF.Next()) {
157             const TopoDS_Shape& aFx=aItLF.Value();
158             if (aFx.IsSame(aF)) {
159               continue;
160             }
161             if (aMCB.Contains(aFx)) {
162               continue;
163             }
164             aMAdd1.Add(aFx);
165           }
166         }//for (; aExp.More(); aExp.Next()){
167         aMCB.Add(aF);
168       }// for (i=1; i<=aNbAdd; ++i) {
169       //
170       aNbAdd1=aMAdd1.Extent();
171       if (!aNbAdd1) {
172         break;// ->make new CB from aMCB
173       }
174       //
175       aMAdd.Clear();
176       for (i=1; i<=aNbAdd1; ++i) {
177         const TopoDS_Shape& aFAdd = aMAdd1(i);
178         aMAdd.Add(aFAdd);
179       }
180     }//while(1) {
181     //
182     aNbF=aMCB.Extent();
183     if (aNbF) {
184       aBB.MakeCompound(aC);
185       //
186       for (i=1; i<=aNbF; ++i) {
187         const TopoDS_Shape& aF=aMCB(i);
188         aBB.Add(aC, aF);  
189         aMP.Add(aF);
190       }
191       theLCB.Append(aC);
192     }
193   }// for (; aIt.More(); aIt.Next()) 
194 }
195 //=======================================================================
196 // function: OrientFacesOnShell
197 // purpose: 
198 //=======================================================================
199   void BOPTools_AlgoTools::OrientFacesOnShell (TopoDS_Shape& aShell)
200 {
201   Standard_Boolean bIsProcessed1, bIsProcessed2;
202   Standard_Integer i, aNbE, aNbF, j;
203   TopAbs_Orientation anOrE1, anOrE2;
204   TopoDS_Face aF1x, aF2x;
205   TopoDS_Shape aShellNew;
206   BOPCol_IndexedDataMapOfShapeListOfShape aEFMap;
207   BOPCol_IndexedMapOfShape aProcessedFaces;
208   BRep_Builder aBB;
209   //
210   BOPTools_AlgoTools::MakeContainer(TopAbs_SHELL, aShellNew);
211   //
212   BOPTools::MapShapesAndAncestors(aShell, TopAbs_EDGE, TopAbs_FACE, aEFMap);
213   aNbE=aEFMap.Extent();
214   // 
215   // One seam edge  in aEFMap contains  2 equivalent faces.
216   for (i=1; i<=aNbE; ++i) {
217     BOPCol_ListOfShape& aLF=aEFMap.ChangeFromIndex(i);
218     aNbF=aLF.Extent();
219     if (aNbF>1) {
220       BOPCol_ListOfShape aLFTmp;
221       BOPCol_IndexedMapOfShape aFM;
222       //
223       BOPCol_ListIteratorOfListOfShape anIt(aLF);
224       for (; anIt.More(); anIt.Next()) {
225         const TopoDS_Shape& aF=anIt.Value();
226         if (!aFM.Contains(aF)) {
227           aFM.Add(aF);
228           aLFTmp.Append(aF);
229         }
230       }
231       aLF.Clear();
232       aLF=aLFTmp;
233     }
234   }
235   //
236   // Do
237   for (i=1; i<=aNbE; ++i) {
238     const TopoDS_Edge& aE=(*(TopoDS_Edge*)(&aEFMap.FindKey(i)));
239     if (BRep_Tool::Degenerated(aE)) {
240       continue;
241     }
242     //
243     const BOPCol_ListOfShape& aLF=aEFMap.FindFromIndex(i);
244     aNbF=aLF.Extent();
245     if (aNbF!=2) {
246       continue;
247     }
248     //
249     TopoDS_Face& aF1=(*(TopoDS_Face*)(&aLF.First()));
250     TopoDS_Face& aF2=(*(TopoDS_Face*)(&aLF.Last()));
251     //    
252     bIsProcessed1=aProcessedFaces.Contains(aF1);
253     bIsProcessed2=aProcessedFaces.Contains(aF2);
254     if (bIsProcessed1 && bIsProcessed2) {
255       continue;
256     }
257     
258     if (!bIsProcessed1 && !bIsProcessed2) {
259       aProcessedFaces.Add(aF1);
260       aBB.Add(aShellNew, aF1);
261       bIsProcessed1=!bIsProcessed1;
262     }
263     //
264     aF1x=aF1;
265     if (bIsProcessed1) {
266       j=aProcessedFaces.FindIndex(aF1);
267       aF1x=(*(TopoDS_Face*)(&aProcessedFaces.FindKey(j)));
268     }
269     //
270     aF2x=aF2;
271     if (bIsProcessed2) {
272       j=aProcessedFaces.FindIndex(aF2);
273       aF2x=(*(TopoDS_Face*)(&aProcessedFaces.FindKey(j)));
274     }
275     //
276     anOrE1=Orientation(aE, aF1x); 
277     anOrE2=Orientation(aE, aF2x);
278     //
279     if (bIsProcessed1 && !bIsProcessed2) {
280       if (anOrE1==anOrE2) {
281         if (!BRep_Tool::IsClosed(aE, aF1) &&
282             !BRep_Tool::IsClosed(aE, aF2)) {
283           aF2.Reverse();
284         }
285       }
286       aProcessedFaces.Add(aF2);
287       aBB.Add(aShellNew, aF2);
288     }
289     else if (!bIsProcessed1 && bIsProcessed2) {
290       if (anOrE1==anOrE2) {
291         if (!BRep_Tool::IsClosed(aE, aF1) &&
292             !BRep_Tool::IsClosed(aE, aF2)) {
293           aF1.Reverse();
294         }
295       }
296       aProcessedFaces.Add(aF1);
297       aBB.Add(aShellNew, aF1);
298     }
299   }
300   //
301   //
302   for (i=1; i<=aNbE; ++i) {
303     const TopoDS_Edge& aE=(*(TopoDS_Edge*)(&aEFMap.FindKey(i)));
304     if (BRep_Tool::Degenerated(aE)) {
305       continue;
306     }
307     //
308     const BOPCol_ListOfShape& aLF=aEFMap.FindFromIndex(i);
309     aNbF=aLF.Extent();
310     if (aNbF!=2) {
311       BOPCol_ListIteratorOfListOfShape anIt(aLF);
312       for(; anIt.More(); anIt.Next()) {
313         const TopoDS_Face& aF=(*(TopoDS_Face*)(&anIt.Value()));
314         if (!aProcessedFaces.Contains(aF)) {
315           aProcessedFaces.Add(aF);
316           aBB.Add(aShellNew, aF);
317         }
318       }
319     }
320   }
321   aShell=aShellNew;
322 }
323 //=======================================================================
324 //function : Orientation
325 //purpose  :
326 //=======================================================================
327   TopAbs_Orientation Orientation(const TopoDS_Edge& anE,
328                                  const TopoDS_Face& aF)
329 {
330   TopAbs_Orientation anOr=TopAbs_INTERNAL;
331
332   TopExp_Explorer anExp;
333   anExp.Init(aF, TopAbs_EDGE);
334   for (; anExp.More(); anExp.Next()) {
335     const TopoDS_Edge& anEF1=(*(TopoDS_Edge*)(&anExp.Current()));
336     if (anEF1.IsSame(anE)) {
337       anOr=anEF1.Orientation();
338       break;
339     }
340   }
341   return anOr;
342 }
343
344 ////////////
345
346
347 //=======================================================================
348 // function: MakeConnexityBlock.
349 // purpose: 
350 //=======================================================================
351   void BOPTools_AlgoTools::MakeConnexityBlock (BOPCol_ListOfShape& theLFIn,
352                                            BOPCol_IndexedMapOfShape& theMEAvoid,
353                                            BOPCol_ListOfShape& theLCB,
354                                            const Handle(NCollection_BaseAllocator)& theAllocator)
355 {
356   Standard_Integer  aNbF, aNbAdd1, aNbAdd, i;
357   TopExp_Explorer aExp;
358   BOPCol_ListIteratorOfListOfShape aIt;
359   //
360   BOPCol_IndexedMapOfShape aMCB(100, theAllocator);
361   BOPCol_IndexedMapOfShape aMAdd(100, theAllocator);
362   BOPCol_IndexedMapOfShape aMAdd1(100, theAllocator);
363   BOPCol_IndexedDataMapOfShapeListOfShape aMEF(100, theAllocator);
364   //
365   // 1. aMEF
366   aNbF=theLFIn.Extent();
367   aIt.Initialize(theLFIn);
368   for (; aIt.More(); aIt.Next()) {
369     const TopoDS_Shape& aF=aIt.Value();      
370     BOPTools::MapShapesAndAncestors(aF, TopAbs_EDGE, TopAbs_FACE, aMEF);
371   }
372   //
373   // 2. aMCB
374   const TopoDS_Shape& aF1=theLFIn.First();
375   aMAdd.Add(aF1);
376   //
377   for(;;) {
378     aMAdd1.Clear();
379     aNbAdd = aMAdd.Extent();
380     for (i=1; i<=aNbAdd; ++i) {
381       const TopoDS_Shape& aF=aMAdd(i);
382       //
383       //aMAdd1.Clear();
384       aExp.Init(aF, TopAbs_EDGE);
385       for (; aExp.More(); aExp.Next()) {
386         const TopoDS_Shape& aE=aExp.Current();
387         if (theMEAvoid.Contains(aE)){
388           continue;
389         }
390         //
391         const BOPCol_ListOfShape& aLF=aMEF.FindFromKey(aE);
392         aIt.Initialize(aLF);
393         for (; aIt.More(); aIt.Next()) {
394           const TopoDS_Shape& aFx=aIt.Value();
395           if (aFx.IsSame(aF)) {
396             continue;
397           }
398           if (aMCB.Contains(aFx)) {
399             continue;
400           }
401           aMAdd1.Add(aFx);
402         }
403       }//for (; aExp.More(); aExp.Next()){
404       aMCB.Add(aF);
405     }// for (i=1; i<=aNbAdd; ++i) {
406     //
407     aNbAdd1=aMAdd1.Extent();
408     if (!aNbAdd1) {
409       break;
410     }
411     //
412     aMAdd.Clear();
413     for (i=1; i<=aNbAdd1; ++i) {
414       const TopoDS_Shape& aFAdd=aMAdd1(i);
415       aMAdd.Add(aFAdd);
416     }
417     //
418   }//while(1) {
419   
420   //
421   aNbF=aMCB.Extent();
422   for (i=1; i<=aNbF; ++i) {
423     const TopoDS_Shape& aF=aMCB(i);
424     theLCB.Append(aF);
425   }
426 }
427 //=======================================================================
428 // function:  ComputeStateByOnePoint
429 // purpose: 
430 //=======================================================================
431   TopAbs_State BOPTools_AlgoTools::ComputeStateByOnePoint(const TopoDS_Shape& theS,
432                                                       const TopoDS_Solid& theRef,
433                                                       const Standard_Real theTol,
434                                                       Handle(BOPInt_Context)& theContext)
435 {
436   TopAbs_State aState;
437   TopAbs_ShapeEnum aType;
438   //
439   aState=TopAbs_UNKNOWN;
440   aType=theS.ShapeType();
441   if (aType==TopAbs_VERTEX) {
442     const TopoDS_Vertex& aV=(*(TopoDS_Vertex*)(&theS));
443     aState=BOPTools_AlgoTools::ComputeState(aV, theRef, theTol, theContext);
444   }
445   else if (aType==TopAbs_EDGE) {
446     const TopoDS_Edge& aE=(*(TopoDS_Edge*)(&theS));
447     aState=BOPTools_AlgoTools::ComputeState(aE, theRef, theTol, theContext);
448   }
449   return aState;
450 }
451
452 //=======================================================================
453 // function:  ComputeState
454 // purpose: 
455 //=======================================================================
456   TopAbs_State BOPTools_AlgoTools::ComputeState(const TopoDS_Face& theF,
457                                             const TopoDS_Solid& theRef,
458                                             const Standard_Real theTol,
459                                             BOPCol_IndexedMapOfShape& theBounds,
460                                             Handle(BOPInt_Context)& theContext)
461 {
462   TopAbs_State aState;
463   TopExp_Explorer aExp; 
464   TopoDS_Edge aE1;
465   gp_Pnt2d aP2D;
466   gp_Pnt aP3D; 
467   //
468   aState=TopAbs_UNKNOWN;
469   //
470   aExp.Init(theF, TopAbs_EDGE);
471   for (; aExp.More(); aExp.Next()) {
472     const TopoDS_Edge& aSE=(*(TopoDS_Edge*)(&aExp.Current()));
473     if (BRep_Tool::Degenerated(aSE)) {
474       continue;
475     }
476     //
477     if (!theBounds.Contains(aSE)) {
478       const TopoDS_Edge& aE=(*(TopoDS_Edge*)(&aSE));
479       aState=BOPTools_AlgoTools::ComputeState(aE, theRef, theTol, theContext);
480       return aState;
481     }
482     if (aE1.IsNull()) {
483       aE1=(*(TopoDS_Edge*)(&aSE));
484     }
485   }
486   // !!<- process edges that are all on theRef
487   if (!aE1.IsNull()) {
488     BOPTools_AlgoTools3D::PointNearEdge(aE1, theF, aP2D, aP3D, theContext);
489     aState=BOPTools_AlgoTools::ComputeState(aP3D, theRef, theTol, theContext);
490   }
491   //
492   return aState;
493 }
494 //=======================================================================
495 // function:  ComputeState
496 // purpose: 
497 //=======================================================================
498   TopAbs_State BOPTools_AlgoTools::ComputeState(const TopoDS_Vertex& theV,
499                                             const TopoDS_Solid& theRef,
500                                             const Standard_Real theTol,
501                                             Handle(BOPInt_Context)& theContext)
502 {
503   TopAbs_State aState;
504   gp_Pnt aP3D; 
505   //
506   aP3D=BRep_Tool::Pnt(theV);
507   aState=BOPTools_AlgoTools::ComputeState(aP3D, theRef, theTol, theContext);
508   return aState;
509 }
510 //=======================================================================
511 // function:  ComputeState
512 // purpose: 
513 //=======================================================================
514   TopAbs_State BOPTools_AlgoTools::ComputeState(const TopoDS_Edge& theE,
515                                             const TopoDS_Solid& theRef,
516                                             const Standard_Real theTol,
517                                             Handle(BOPInt_Context)& theContext)
518 {
519   Standard_Real aT1, aT2, aT = 0.;
520   TopAbs_State aState;
521   Handle(Geom_Curve) aC3D;
522   gp_Pnt aP3D; 
523   //
524   aC3D = BRep_Tool::Curve(theE, aT1, aT2);
525   //
526   if(aC3D.IsNull()) {
527     //it means that we are in degenerated edge
528     const TopoDS_Vertex& aV = TopExp::FirstVertex(theE);
529     if(aV.IsNull()){
530       return TopAbs_UNKNOWN;
531     }
532     aP3D=BRep_Tool::Pnt(aV);
533   }
534   else {//usual case
535     Standard_Boolean bF2Inf, bL2Inf;
536     Standard_Real dT=10.;
537     //
538     bF2Inf = Precision::IsNegativeInfinite(aT1);
539     bL2Inf = Precision::IsPositiveInfinite(aT2);
540     //
541     if (bF2Inf && !bL2Inf) {
542       aT=aT2-dT;
543     }
544     else if (!bF2Inf && bL2Inf) {
545       aT=aT1+dT;
546     }
547     else if (bF2Inf && bL2Inf) {
548       aT=0.;
549     }
550     else {
551       aT=IntTools_Tools::IntermediatePoint(aT1, aT2);
552     }
553     aC3D->D0(aT, aP3D);
554   }
555   //
556   aState=BOPTools_AlgoTools::ComputeState(aP3D, theRef, theTol, theContext);
557   //
558   return aState;
559 }
560 //=======================================================================
561 // function:  ComputeState
562 // purpose: 
563 //=======================================================================
564   TopAbs_State BOPTools_AlgoTools::ComputeState(const gp_Pnt& theP,
565                                             const TopoDS_Solid& theRef,
566                                             const Standard_Real theTol,
567                                             Handle(BOPInt_Context)& theContext)
568 {
569   TopAbs_State aState;
570   //
571   BRepClass3d_SolidClassifier& aSC=theContext->SolidClassifier(theRef);
572   aSC.Perform(theP, theTol);
573   //
574   aState=aSC.State();
575   //
576   return aState;
577 }
578 //=======================================================================
579 //function : IsInternalFace
580 //purpose  : 
581 //=======================================================================
582   Standard_Boolean BOPTools_AlgoTools::IsInternalFace(const TopoDS_Face& theFace,
583                                                       const TopoDS_Solid& theSolid,
584                                                       BOPCol_IndexedDataMapOfShapeListOfShape& theMEF,
585                                                       const Standard_Real theTol,
586                                                       Handle(BOPInt_Context)& theContext)
587 {
588   Standard_Boolean bRet, bDegenerated;
589   Standard_Integer aNbF;
590   TopAbs_Orientation aOr;
591   TopoDS_Edge aEL;
592   TopExp_Explorer aExp;
593   BOPCol_ListIteratorOfListOfShape aItF;
594   //
595   bRet=Standard_False;
596   //
597   // 1 Try to find an edge from theFace in theMEF
598   aExp.Init(theFace, TopAbs_EDGE);
599   for(; aExp.More(); aExp.Next()) {
600     const TopoDS_Edge& aE=(*(TopoDS_Edge*)(&aExp.Current()));
601     if (!theMEF.Contains(aE)) {
602       continue;
603     }
604     //
605     aOr=aE.Orientation();
606     if (aOr==TopAbs_INTERNAL) {
607       continue;
608     }
609     bDegenerated=BRep_Tool::Degenerated(aE);
610     if (bDegenerated){
611       continue;
612     }
613     // aE
614     BOPCol_ListOfShape& aLF=theMEF.ChangeFromKey(aE);
615     aNbF=aLF.Extent();
616     if (!aNbF) {
617       return bRet; // it can not be so
618     }
619     else if (aNbF==1) {
620       // aE is internal edge on aLF.First()
621       const TopoDS_Face& aF1=(*(TopoDS_Face*)(&aLF.First()));
622       bRet=BOPTools_AlgoTools::IsInternalFace(theFace, aE, aF1, aF1, theContext);
623       return bRet;
624     }
625     else if (aNbF==2) {
626       const TopoDS_Face& aF1=(*(TopoDS_Face*)(&aLF.First()));
627       const TopoDS_Face& aF2=(*(TopoDS_Face*)(&aLF.Last()));
628       //
629       if (aF2.IsSame(aF1) && BRep_Tool::IsClosed(aE, aF1)) {
630         // treat as it was for 1 face
631         bRet=BOPTools_AlgoTools::IsInternalFace(theFace, aE, aF1, aF2, theContext);
632         return bRet;
633       }
634     }
635     if (aNbF%2) {
636       return bRet; // it can not be so
637     }
638     else { // aNbF=2,4,6,8,...
639       bRet=BOPTools_AlgoTools::IsInternalFace(theFace, aE, aLF, theContext);
640       return bRet;
641     }
642   }//for(; aExp.More(); aExp.Next()) {
643   //
644   //========================================
645   // 2. Classify face using classifier
646   //
647   TopAbs_State aState;
648   BOPCol_IndexedMapOfShape aBounds;
649   //
650   aState=BOPTools_AlgoTools::ComputeState(theFace, theSolid, theTol, aBounds, theContext);
651   bRet=(aState==TopAbs_IN);
652   //
653   return bRet;
654 }
655 //=======================================================================
656 //function : IsInternalFace
657 //purpose  : 
658 //=======================================================================
659   Standard_Boolean BOPTools_AlgoTools::IsInternalFace(const TopoDS_Face& theFace,
660                                                       const TopoDS_Edge& theEdge,
661                                                       BOPCol_ListOfShape& theLF,
662                                                       Handle(BOPInt_Context)& theContext)
663 {
664   Standard_Boolean bRet;
665   Standard_Boolean aNbF;
666   //
667   bRet=Standard_False;
668   //
669   aNbF=theLF.Extent();
670   if (aNbF==2) {
671     const TopoDS_Face& aF1=(*(TopoDS_Face*)(&theLF.First()));
672     const TopoDS_Face& aF2=(*(TopoDS_Face*)(&theLF.Last()));
673     bRet=BOPTools_AlgoTools::IsInternalFace(theFace, theEdge, aF1, aF2, theContext);
674     return bRet;
675   }
676   //
677   else {
678     BOPTools_ListOfCoupleOfShape aLCFF;
679     BOPTools_ListIteratorOfListOfCoupleOfShape aIt;
680     //
681     FindFacePairs(theEdge, theLF, aLCFF, theContext);
682     //
683     aIt.Initialize(aLCFF);
684     for (; aIt.More(); aIt.Next()) {
685       BOPTools_CoupleOfShape& aCSFF=aIt.ChangeValue();
686       //
687       const TopoDS_Face& aF1=(*(TopoDS_Face*)(&aCSFF.Shape1()));
688       const TopoDS_Face& aF2=(*(TopoDS_Face*)(&aCSFF.Shape2()));
689       bRet=BOPTools_AlgoTools::IsInternalFace(theFace, theEdge, aF1, aF2, theContext);
690       if (bRet) {
691         return bRet;
692       }
693     }
694   }
695   return bRet;
696 }
697 //=======================================================================
698 //function : IsInternalFace
699 //purpose  : 
700 //=======================================================================
701   Standard_Boolean BOPTools_AlgoTools::IsInternalFace(const TopoDS_Face& theFace,
702                                                       const TopoDS_Edge& theEdge,
703                                                       const TopoDS_Face& theFace1,
704                                                       const TopoDS_Face& theFace2,
705                                                       Handle(BOPInt_Context)& theContext)
706 {
707   Standard_Boolean bRet;
708   TopoDS_Edge aE1, aE2;
709   TopoDS_Face aFOff;
710   BOPTools_ListOfCoupleOfShape theLCSOff;
711   BOPTools_CoupleOfShape aCS1, aCS2;
712   //
713   BOPTools_AlgoTools::GetEdgeOnFace(theEdge, theFace1, aE1);
714   if (aE1.Orientation()==TopAbs_INTERNAL) {
715     aE2=aE1;
716     aE1.Orientation(TopAbs_FORWARD);
717     aE2.Orientation(TopAbs_REVERSED);
718   }
719   else if (theFace1==theFace2) {
720     aE2=aE1;
721     aE1.Orientation(TopAbs_FORWARD);
722     aE2.Orientation(TopAbs_REVERSED);
723   }
724   else {
725     BOPTools_AlgoTools::GetEdgeOnFace(theEdge, theFace2, aE2);
726   }
727   //
728   aCS1.SetShape1(theEdge);
729   aCS1.SetShape2(theFace);
730   theLCSOff.Append(aCS1);
731   //
732   aCS2.SetShape1(aE2);
733   aCS2.SetShape2(theFace2);
734   theLCSOff.Append(aCS2);
735   //
736   GetFaceOff(aE1, theFace1, theLCSOff, aFOff, theContext);
737   //
738   bRet = theFace.IsEqual(aFOff);
739   return bRet;
740 }
741 //=======================================================================
742 //function : GetFaceOff
743 //purpose  : 
744 //=======================================================================
745   void BOPTools_AlgoTools::GetFaceOff(const TopoDS_Edge& theE1,
746                                       const TopoDS_Face& theF1,
747                                       BOPTools_ListOfCoupleOfShape& theLCSOff,
748                                       TopoDS_Face& theFOff,
749                                       Handle(BOPInt_Context)& theContext)
750 {
751   Standard_Real aT, aT1, aT2, aAngle, aTwoPI, aAngleMin;
752   gp_Pnt aPn1, aPn2, aPx;
753   gp_Dir aDN1, aDN2, aDBF, aDBF2, aDTF;
754   gp_Vec aVTgt;
755   TopAbs_Orientation aOr;
756   Handle(Geom_Curve)aC3D;
757   BOPTools_ListIteratorOfListOfCoupleOfShape aIt;
758   //
759   aAngleMin=100.;
760   aTwoPI=M_PI+M_PI;
761   aC3D =BRep_Tool::Curve(theE1, aT1, aT2);
762   aT=BOPTools_AlgoTools2D::IntermediatePoint(aT1, aT2);
763   aC3D->D0(aT, aPx);
764   //
765   BOPTools_AlgoTools2D::EdgeTangent(theE1, aT, aVTgt);
766   gp_Dir aDTgt(aVTgt), aDTgt2;
767   aOr = theE1.Orientation();
768   //
769   GetFaceDir(theE1, theF1, aPx, aT, aDTgt, aDN1, aDBF, theContext);
770   //
771   aDTF=aDN1^aDBF;
772   //
773   aIt.Initialize(theLCSOff);
774   for (; aIt.More(); aIt.Next()) {
775     const BOPTools_CoupleOfShape& aCS=aIt.Value();
776     const TopoDS_Edge& aE2=(*(TopoDS_Edge*)(&aCS.Shape1()));
777     const TopoDS_Face& aF2=(*(TopoDS_Face*)(&aCS.Shape2()));
778     //
779     aDTgt2 = (aE2.Orientation()==aOr) ? aDTgt : aDTgt.Reversed();
780     GetFaceDir(aE2, aF2, aPx, aT, aDTgt2, aDN2, aDBF2, theContext);
781     //Angle
782     aAngle=AngleWithRef(aDBF, aDBF2, aDTF);
783     //
784     if(aAngle<0.) {
785       aAngle=aTwoPI+aAngle;
786     }
787     //
788     if (aAngle<Precision::Angular()) {
789       if (aF2==theF1) {
790         aAngle=M_PI;
791       }
792       else if (aF2.IsSame(theF1)) {
793         aAngle=aTwoPI;
794       }
795     }
796     //
797     if (aAngle<aAngleMin){
798       aAngleMin=aAngle;
799       theFOff=aF2;
800     }
801   }
802 }
803 //=======================================================================
804 //function : GetEdgeOff
805 //purpose  : 
806 //=======================================================================
807   Standard_Boolean BOPTools_AlgoTools::GetEdgeOff(const TopoDS_Edge& theE1,
808                                               const TopoDS_Face& theF2,
809                                               TopoDS_Edge& theE2)
810 {
811   Standard_Boolean bFound;
812   TopAbs_Orientation aOr1, aOr1C, aOr2;
813   TopExp_Explorer anExp;
814   //
815   bFound=Standard_False;
816   aOr1=theE1.Orientation();
817   aOr1C=TopAbs::Reverse(aOr1);
818   //
819   anExp.Init(theF2, TopAbs_EDGE);
820   for (; anExp.More(); anExp.Next()) {
821     const TopoDS_Edge& aEF2=(*(TopoDS_Edge*)(&anExp.Current()));
822     if (aEF2.IsSame(theE1)) {
823       aOr2=aEF2.Orientation();
824       if (aOr2==aOr1C) {
825         theE2=aEF2;
826         bFound=!bFound;
827         return bFound;
828       }
829     }
830   }
831   return bFound;
832 }
833
834 //=======================================================================
835 //function : AreFacesSameDomain
836 //purpose  : 
837 //=======================================================================
838   Standard_Boolean BOPTools_AlgoTools::AreFacesSameDomain(const TopoDS_Face& theF1,
839                                                           const TopoDS_Face& theF2,
840                                                           Handle(BOPInt_Context)& theContext)
841 {
842   Standard_Boolean bFlag;
843   Standard_Integer iErr;
844   Standard_Real aTolF1, aTolF2, aTol;
845   gp_Pnt2d aP2D;
846   gp_Pnt aP;
847   TopoDS_Face aF1, aF2;
848   TopoDS_Edge aE1;
849   TopExp_Explorer aExp;
850   //
851   bFlag=Standard_False;
852   //
853   aF1=theF1;
854   aF1.Orientation(TopAbs_FORWARD);
855   aF2=theF2;
856   aF2.Orientation(TopAbs_FORWARD);
857   //
858   aTolF1=BRep_Tool::Tolerance(aF1);
859   // 1
860   aExp.Init(aF1, TopAbs_EDGE);
861   for (; aExp.More(); aExp.Next()) {
862     aE1=(*(TopoDS_Edge*)(&aExp.Current()));
863     if (!BRep_Tool::Degenerated(aE1)) {
864       Standard_Real aTolE = BRep_Tool::Tolerance(aE1);
865       aTolF1 = (aTolE > aTolF1) ? aTolE : aTolF1;
866     }
867   }
868   // 2
869   aTolF2=BRep_Tool::Tolerance(aF2);
870   aTol=aTolF1+aTolF2;
871   //
872   iErr = BOPTools_AlgoTools3D::PointInFace(aF1, aP, aP2D, theContext);
873   if (!iErr) {
874     bFlag=theContext->IsValidPointForFace(aP, aF2, aTol);
875   }
876   //
877   return bFlag;
878 }
879
880 //=======================================================================
881 //function : CheckSameGeom
882 //purpose  : 
883 //=======================================================================
884   Standard_Boolean BOPTools_AlgoTools::CheckSameGeom(const TopoDS_Face& theF1,
885                                                  const TopoDS_Face& theF2,
886                                                  Handle(BOPInt_Context)& theContext)
887 {
888   Standard_Boolean bRet;
889   Standard_Real aTolF1, aTolF2, aTol;
890   gp_Pnt2d aP2D;
891   gp_Pnt aP;
892   TopExp_Explorer aExp;
893   //
894   bRet=Standard_False;
895   aExp.Init(theF1, TopAbs_EDGE);
896   for (; aExp.More(); aExp.Next()) {
897     const TopoDS_Edge& aE=(*(TopoDS_Edge*)(&aExp.Current()));
898     if (!BRep_Tool::Degenerated(aE)) {
899       aTolF1=BRep_Tool::Tolerance(theF1);
900       aTolF2=BRep_Tool::Tolerance(theF2);
901       aTol=aTolF1+aTolF2;
902       BOPTools_AlgoTools3D::PointNearEdge(aE, theF1, aP2D, aP, theContext);
903       bRet=theContext->IsValidPointForFace(aP, theF2, aTol);
904       break;
905     }
906   }
907   return bRet;
908 }
909 //=======================================================================
910 // function: Sense
911 // purpose: 
912 //=======================================================================
913   Standard_Integer BOPTools_AlgoTools::Sense (const TopoDS_Face& theF1,
914                                           const TopoDS_Face& theF2)
915 {
916   Standard_Integer iSense=0;
917   gp_Dir aDNF1, aDNF2;
918   TopoDS_Edge aE1, aE2;
919   TopExp_Explorer aExp;
920   //
921   aExp.Init(theF1, TopAbs_EDGE);
922   for (; aExp.More(); aExp.Next()) {
923     aE1=(*(TopoDS_Edge*)(&aExp.Current()));
924     if (!BRep_Tool::Degenerated(aE1)) {
925       if (!BRep_Tool::IsClosed(aE1, theF1)) {
926         break;
927       }
928     }
929   }
930   //
931   aExp.Init(theF2, TopAbs_EDGE);
932   for (; aExp.More(); aExp.Next()) {
933     aE2=(*(TopoDS_Edge*)(&aExp.Current()));
934     if (!BRep_Tool::Degenerated(aE2)) {
935       if (!BRep_Tool::IsClosed(aE2, theF2)) {
936         if (aE2.IsSame(aE1)) {
937           iSense=1;
938           break;
939         }
940       }
941     }
942   }
943   //
944   if (!iSense) {
945     return iSense;
946   }
947   //
948   BOPTools_AlgoTools3D::GetNormalToFaceOnEdge(aE1, theF1, aDNF1);
949   BOPTools_AlgoTools3D::GetNormalToFaceOnEdge(aE2, theF2, aDNF2);
950   //
951   iSense=BOPTools_AlgoTools3D::SenseFlag(aDNF1, aDNF2);
952   //
953   return iSense;
954 }
955 //=======================================================================
956 // function: IsSplitToReverse
957 // purpose: 
958 //=======================================================================
959   Standard_Boolean BOPTools_AlgoTools::IsSplitToReverse(const TopoDS_Shape& theSp,
960                                                     const TopoDS_Shape& theSr,
961                                                     Handle(BOPInt_Context)& theContext)
962 {
963   Standard_Boolean bRet;
964   TopAbs_ShapeEnum aType;
965   //
966   bRet=Standard_False;
967   //
968   aType=theSp.ShapeType();
969   switch (aType) {
970     case TopAbs_EDGE: {
971       const TopoDS_Edge& aESp=(*(TopoDS_Edge*)(&theSp));
972       const TopoDS_Edge& aESr=(*(TopoDS_Edge*)(&theSr));
973       bRet=BOPTools_AlgoTools::IsSplitToReverse(aESp, aESr, theContext);
974     }
975       break;
976       //
977     case TopAbs_FACE: {
978       const TopoDS_Face& aFSp=(*(TopoDS_Face*)(&theSp));
979       const TopoDS_Face& aFSr=(*(TopoDS_Face*)(&theSr));
980       bRet=BOPTools_AlgoTools::IsSplitToReverse(aFSp, aFSr, theContext);
981     }
982       break;
983       //
984     default:
985       break;
986   }
987   return bRet;
988 }
989 //=======================================================================
990 //function :IsSplitToReverse
991 //purpose  : 
992 //=======================================================================
993   Standard_Boolean BOPTools_AlgoTools::IsSplitToReverse(const TopoDS_Face& theFSp,
994                                                     const TopoDS_Face& theFSr,
995                                                     Handle(BOPInt_Context)& theContext)
996 {
997   Standard_Boolean bRet, bFound, bInFace;
998   Standard_Real aT1, aT2, aT, aU, aV, aScPr;
999   gp_Pnt aPFSp, aPFSr;
1000   gp_Dir aDNFSp;
1001   gp_Vec aD1U, aD1V;
1002   Handle(Geom_Surface) aSr, aSp;
1003   TopAbs_Orientation aOrSr, aOrSp;
1004   TopExp_Explorer anExp;
1005   TopoDS_Edge aESp;
1006   //
1007   bRet=Standard_False;
1008   //
1009   aSr=BRep_Tool::Surface(theFSr);
1010   aSp=BRep_Tool::Surface(theFSp);
1011   if (aSr==aSp) {
1012     aOrSr=theFSr.Orientation();
1013     aOrSp=theFSp.Orientation();
1014     bRet=(aOrSr!=aOrSp);
1015     return bRet;
1016   }
1017   //
1018   bFound=Standard_False;
1019   anExp.Init(theFSp, TopAbs_EDGE);
1020   for (; anExp.More(); anExp.Next()) {
1021     aESp=(*(TopoDS_Edge*)(&anExp.Current()));
1022     if (!BRep_Tool::Degenerated(aESp)) {
1023       if (!BRep_Tool::IsClosed(aESp, theFSp)) {
1024         bFound=!bFound;
1025         break;
1026       }
1027     }
1028   }
1029   if (!bFound) {
1030     Standard_Boolean bFlag;
1031     Standard_Integer iErr;
1032     gp_Pnt2d aP2DFSp;
1033     //
1034     iErr=BOPTools_AlgoTools3D::PointInFace(theFSp, aPFSp, aP2DFSp, theContext);
1035     if (iErr) {
1036       return bRet;
1037     }
1038     //
1039     aP2DFSp.Coord(aU, aV);
1040     bFlag=BOPTools_AlgoTools3D::GetNormalToSurface(aSp, aU, aV, aDNFSp);
1041     if (!bFlag) {
1042       return bRet;
1043     }
1044   }
1045   else {
1046     BRep_Tool::Range(aESp, aT1, aT2);
1047     aT=BOPTools_AlgoTools2D::IntermediatePoint(aT1, aT2);
1048     BOPTools_AlgoTools3D::GetApproxNormalToFaceOnEdge(aESp, theFSp, aT, aPFSp, aDNFSp, theContext);
1049   }
1050   //
1051   // Parts of theContext->ComputeVS(..) 
1052   GeomAPI_ProjectPointOnSurf& aProjector=theContext->ProjPS(theFSr);
1053   aProjector.Perform(aPFSp);
1054   if (!aProjector.IsDone()) {
1055     return bRet;
1056   }
1057   //
1058   aProjector.LowerDistanceParameters(aU, aV);
1059   gp_Pnt2d aP2D(aU, aV);
1060   bInFace=theContext->IsPointInFace (theFSr, aP2D);
1061   if (!bInFace) {
1062     return bRet;
1063   }
1064   //
1065   aSr->D1(aU, aV, aPFSr, aD1U, aD1V);
1066   gp_Dir aDD1U(aD1U); 
1067   gp_Dir aDD1V(aD1V);
1068   gp_Dir aDNFSr=aDD1U^aDD1V; 
1069   if (theFSr.Orientation()==TopAbs_REVERSED){
1070     aDNFSr.Reverse();
1071   }
1072   //
1073   aScPr=aDNFSp*aDNFSr;
1074   bRet=(aScPr<0.);
1075   //
1076   return bRet;
1077 }
1078 //=======================================================================
1079 //function :IsSplitToReverse
1080 //purpose  : 
1081 //=======================================================================
1082   Standard_Boolean BOPTools_AlgoTools::IsSplitToReverse(const TopoDS_Edge& aEF1,
1083                                                     const TopoDS_Edge& aEF2,
1084                                                     Handle(BOPInt_Context)& theContext)
1085 {
1086   Standard_Boolean bRet, bIsDegenerated;
1087   //
1088   bRet=Standard_False;
1089   bIsDegenerated=(BRep_Tool::Degenerated(aEF1) || 
1090                   BRep_Tool::Degenerated(aEF2));
1091   if (bIsDegenerated) {
1092     return bRet;
1093   }
1094   //
1095   Standard_Real a, b;
1096   TopAbs_Orientation aOrE, aOrSp;
1097   Handle(Geom_Curve)aC1, aC2;
1098   //
1099   aC2=BRep_Tool::Curve(aEF2, a, b);
1100   aC1=BRep_Tool::Curve(aEF1, a, b);
1101   //
1102   if (aC1==aC2) {
1103     aOrE=aEF2.Orientation();
1104     aOrSp=aEF1.Orientation();
1105     bRet=(aOrE!=aOrSp);
1106     return bRet;
1107   }
1108   //
1109   Standard_Real aT1, aT2, aScPr;
1110   gp_Vec aV1, aV2;
1111   gp_Pnt aP;
1112   //
1113   aT1=BOPTools_AlgoTools2D::IntermediatePoint(a, b);
1114   aC1->D0(aT1, aP);
1115   BOPTools_AlgoTools2D::EdgeTangent(aEF1, aT1, aV1);
1116   gp_Dir aDT1(aV1);
1117   //
1118   theContext->ProjectPointOnEdge(aP, aEF2, aT2);
1119   //
1120   BOPTools_AlgoTools2D::EdgeTangent(aEF2, aT2, aV2);
1121   gp_Dir aDT2(aV2);
1122   //
1123   aScPr=aDT1*aDT2;
1124   bRet=(aScPr<0.);
1125   //
1126   return bRet;
1127 }
1128
1129 //=======================================================================
1130 //function : IsHole
1131 //purpose  : 
1132 //=======================================================================
1133   Standard_Boolean BOPTools_AlgoTools::IsHole(const TopoDS_Shape& aW,
1134                                           const TopoDS_Shape& aFace)
1135 {
1136   Standard_Boolean bIsHole;
1137   Standard_Integer i, aNbS;
1138   Standard_Real aT1, aT2, aS;
1139   Standard_Real aU1, aU2, aU, dU;
1140   Standard_Real aX1, aY1, aX0, aY0;
1141   TopAbs_Orientation aOr;
1142   
1143   gp_Pnt2d aP2D0, aP2D1;
1144   Handle(Geom2d_Curve) aC2D;
1145   TopoDS_Face aF, aFF;
1146   TopoDS_Iterator aItW;
1147   //
1148   bIsHole=Standard_False;
1149   //
1150   aF=(*(TopoDS_Face *)(&aFace));
1151   aFF=aF;
1152   aFF.Orientation(TopAbs_FORWARD);
1153   //
1154   aS=0.;
1155   aItW.Initialize(aW);
1156   for (; aItW.More(); aItW.Next()) { 
1157     const TopoDS_Edge& aE=(*(TopoDS_Edge *)(&aItW.Value()));
1158     aOr=aE.Orientation();
1159     if (!(aOr==TopAbs_FORWARD || 
1160           aOr==TopAbs_REVERSED)) {
1161       continue;
1162     }
1163     //
1164     aC2D=BRep_Tool::CurveOnSurface(aE, aFF, aT1, aT2);
1165     if (aC2D.IsNull()) {
1166       break; //xx
1167     }
1168     //
1169     BRepAdaptor_Curve2d aBAC2D(aE, aFF);
1170     aNbS=Geom2dInt_Geom2dCurveTool::NbSamples(aBAC2D);
1171     if (aNbS>2) {
1172       aNbS*=4;
1173     }
1174     //
1175     dU=(aT2-aT1)/(Standard_Real)(aNbS-1);
1176     aU =aT1;
1177     aU1=aT1;
1178     aU2=aT2;
1179     if (aOr==TopAbs_REVERSED) {
1180       aU =aT2;
1181       aU1=aT2;
1182       aU2=aT1;
1183       dU=-dU;
1184     }
1185     //
1186     aC2D->D0(aU, aP2D0);
1187     for(i=2; i<=aNbS; i++) {
1188       aU=aU1+(i-1)*dU;
1189       aC2D->D0(aU, aP2D1);
1190       aP2D0.Coord(aX0, aY0);
1191       aP2D1.Coord(aX1, aY1);
1192       //
1193       aS=aS+(aY0+aY1)*(aX1-aX0); 
1194       //
1195       aP2D0=aP2D1;
1196     }
1197   }//for (; aItW.More(); aItW.Next()) { 
1198   bIsHole=(aS>0.);
1199   return bIsHole;
1200 }
1201
1202 //=======================================================================
1203 // function: MakeContainer
1204 // purpose: 
1205 //=======================================================================
1206   void BOPTools_AlgoTools::MakeContainer(const TopAbs_ShapeEnum theType,
1207                                      TopoDS_Shape& theC)
1208 {
1209   BRep_Builder aBB;
1210   //
1211   switch(theType) {
1212     case TopAbs_COMPOUND:{
1213       TopoDS_Compound aC;
1214       aBB.MakeCompound(aC);
1215       theC=aC;
1216     }
1217       break;
1218       //
1219     case TopAbs_COMPSOLID:{
1220       TopoDS_CompSolid aCS;
1221       aBB.MakeCompSolid(aCS);
1222       theC=aCS;
1223     }
1224       break;
1225       //
1226     case TopAbs_SOLID:{
1227       TopoDS_Solid aSolid;
1228       aBB.MakeSolid(aSolid);
1229       theC=aSolid;
1230     }  
1231       break;
1232       //
1233       //
1234     case TopAbs_SHELL:{
1235       TopoDS_Shell aShell;
1236       aBB.MakeShell(aShell);
1237       theC=aShell;
1238     }  
1239       break;
1240       //
1241     case TopAbs_WIRE: {
1242       TopoDS_Wire aWire;
1243       aBB.MakeWire(aWire);
1244       theC=aWire;
1245     }
1246       break;
1247       //
1248     default:
1249       break;
1250   }
1251 }
1252 //=======================================================================
1253 // function: MakePCurve
1254 // purpose: 
1255 //=======================================================================
1256   void  BOPTools_AlgoTools::MakePCurve(const TopoDS_Edge& aE,
1257                                        const TopoDS_Face& aF1,
1258                                        const TopoDS_Face& aF2,
1259                                        const IntTools_Curve& aIC,
1260                                        const Standard_Boolean bPC1,
1261                                        const Standard_Boolean bPC2)
1262
1263 {
1264   Standard_Integer i;
1265   Standard_Real aTolE, aT1, aT2, aOutFirst, aOutLast, aOutTol;
1266   Handle(Geom2d_Curve) aC2D, aC2DA, aC2Dx1;
1267   TopoDS_Face aFFWD; 
1268   BRep_Builder aBB;
1269   Standard_Boolean bPC;
1270   //
1271   aTolE=BRep_Tool::Tolerance(aE);
1272   //
1273   const Handle(Geom_Curve)& aC3DE=BRep_Tool::Curve(aE, aT1, aT2);
1274   Handle(Geom_TrimmedCurve)aC3DETrim=new Geom_TrimmedCurve(aC3DE, aT1, aT2);
1275   //
1276   for (i=0; i<2; ++i) {
1277     bPC = !i ? bPC1 : bPC2;
1278     if (!bPC) {
1279       continue;
1280     }
1281     //
1282     if (!i) {
1283       aFFWD=aF1;
1284       aC2Dx1=aIC.FirstCurve2d();
1285     }
1286     else {
1287       aFFWD=aF2;
1288       aC2Dx1=aIC.SecondCurve2d();
1289     }
1290     //
1291     aFFWD.Orientation(TopAbs_FORWARD);
1292     //
1293     aC2D=aC2Dx1;
1294     if (aC2D.IsNull()) { 
1295       BOPTools_AlgoTools2D::BuildPCurveForEdgeOnFace(aE, aFFWD);
1296       BOPTools_AlgoTools2D::CurveOnSurface(aE, aFFWD, aC2D, 
1297                                        aOutFirst, aOutLast, 
1298                                        aOutTol);
1299       }
1300     //
1301     if (aC3DE->IsPeriodic()) {
1302       BOPTools_AlgoTools2D::AdjustPCurveOnFace(aFFWD, aT1, aT2,  aC2D, aC2DA); 
1303     }
1304     else {
1305       BOPTools_AlgoTools2D::AdjustPCurveOnFace(aFFWD, aC3DETrim, aC2D, aC2DA); 
1306     }
1307     //
1308     aBB.UpdateEdge(aE, aC2DA, aFFWD, aTolE);
1309     //BRepLib::SameParameter(aE);
1310   }
1311   BRepLib::SameParameter(aE);
1312 }
1313 //=======================================================================
1314 // function: MakeEdge
1315 // purpose: 
1316 //=======================================================================
1317   void  BOPTools_AlgoTools::MakeEdge(const IntTools_Curve& theIC,
1318                                  const TopoDS_Vertex& theV1,
1319                                  const Standard_Real theT1,
1320                                  const TopoDS_Vertex& theV2,
1321                                  const Standard_Real theT2,
1322                                  const Standard_Real theTolR3D,
1323                                  TopoDS_Edge& theE)
1324 {
1325   Standard_Real aTolV;
1326   BRep_Builder aBB;
1327   //
1328   BOPTools_AlgoTools::MakeSectEdge (theIC, theV1, theT1, theV2, theT2, theE);
1329   //
1330   aBB.UpdateEdge(theE, theTolR3D);
1331   //
1332   aTolV=BRep_Tool::Tolerance(theV1);
1333   if (aTolV<theTolR3D) {
1334     aBB.UpdateVertex(theV1, theTolR3D);
1335   }
1336   //
1337   aTolV=BRep_Tool::Tolerance(theV2);
1338   if (aTolV<theTolR3D) {
1339     aBB.UpdateVertex(theV2, theTolR3D);
1340   }
1341 }
1342 //=======================================================================
1343 // function: ComputeVV
1344 // purpose: 
1345 //=======================================================================
1346   Standard_Integer BOPTools_AlgoTools::ComputeVV(const TopoDS_Vertex& aV1, 
1347                                              const gp_Pnt& aP2,
1348                                              const Standard_Real aTolP2)
1349 {
1350   Standard_Real aTolV1, aTolSum, aTolSum2, aD2;
1351   gp_Pnt aP1;
1352   //
1353   aTolV1=BRep_Tool::Tolerance(aV1);
1354   
1355   aTolSum=aTolV1+aTolP2;
1356   aTolSum2=aTolSum*aTolSum;
1357   //
1358   aP1=BRep_Tool::Pnt(aV1);
1359   //
1360   aD2=aP1.SquareDistance(aP2);
1361   if (aD2>aTolSum2) {
1362     return 1;
1363   }
1364   return 0;
1365
1366 //=======================================================================
1367 // function: ComputeVV
1368 // purpose: 
1369 //=======================================================================
1370   Standard_Integer BOPTools_AlgoTools::ComputeVV(const TopoDS_Vertex& aV1, 
1371                                              const TopoDS_Vertex& aV2)
1372 {
1373   Standard_Real aTolV1, aTolV2, aTolSum, aTolSum2, aD2;
1374   gp_Pnt aP1, aP2;
1375   //
1376   aTolV1=BRep_Tool::Tolerance(aV1);
1377   aTolV2=BRep_Tool::Tolerance(aV2);
1378   aTolSum=aTolV1+aTolV2;
1379   aTolSum2=aTolSum*aTolSum;
1380   //
1381   aP1=BRep_Tool::Pnt(aV1);
1382   aP2=BRep_Tool::Pnt(aV2);
1383   //
1384   aD2=aP1.SquareDistance(aP2);
1385   if (aD2>aTolSum2) {
1386     return 1;
1387   }
1388   return 0;
1389 }
1390 //=======================================================================
1391 // function: MakeVertex
1392 // purpose : 
1393 //=======================================================================
1394   void BOPTools_AlgoTools::MakeVertex(BOPCol_ListOfShape& aLV,
1395                                   TopoDS_Vertex& aVnew)
1396 {
1397   Standard_Integer aNb;
1398   Standard_Real aTi, aDi, aDmax;
1399   gp_Pnt aPi, aP;
1400   gp_XYZ aXYZ(0.,0.,0.), aXYZi;
1401   BOPCol_ListIteratorOfListOfShape aIt;
1402   //
1403   aNb=aLV.Extent();
1404   if (aNb) {
1405     aIt.Initialize(aLV);
1406     for (; aIt.More(); aIt.Next()) {
1407       TopoDS_Vertex& aVi=*((TopoDS_Vertex*)(&aIt.Value()));
1408       aPi=BRep_Tool::Pnt(aVi);
1409       aXYZi=aPi.XYZ();
1410       aXYZ=aXYZ+aXYZi;
1411     }
1412     //
1413     aXYZ.Divide((Standard_Real)aNb);
1414     aP.SetXYZ(aXYZ);
1415     //
1416     aDmax=-1.;
1417     aIt.Initialize(aLV);
1418     for (; aIt.More(); aIt.Next()) {
1419       TopoDS_Vertex& aVi=*((TopoDS_Vertex*)(&aIt.Value()));
1420       aPi=BRep_Tool::Pnt(aVi);
1421       aTi=BRep_Tool::Tolerance(aVi);
1422       aDi=aP.SquareDistance(aPi);
1423       aDi=fsqrt(aDi);
1424       aDi=aDi+aTi;
1425       if (aDi > aDmax) {
1426         aDmax=aDi;
1427       }
1428     }
1429     //
1430     BRep_Builder aBB;
1431     aBB.MakeVertex (aVnew, aP, aDmax);
1432   }
1433 }
1434 //=======================================================================
1435 //function : GetEdgeOnFace
1436 //purpose  : 
1437 //=======================================================================
1438   Standard_Boolean BOPTools_AlgoTools::GetEdgeOnFace(const TopoDS_Edge& theE1,
1439                                                  const TopoDS_Face& theF2,
1440                                                  TopoDS_Edge& theE2)
1441 {
1442   Standard_Boolean bFound;
1443   TopoDS_Iterator aItF, aItW;
1444   //
1445   bFound=Standard_False;
1446   //
1447   aItF.Initialize(theF2);
1448   for (; aItF.More(); aItF.Next()) {
1449     const TopoDS_Shape& aW=aItF.Value();
1450     aItW.Initialize(aW);
1451     for (; aItW.More(); aItW.Next()) {
1452       const TopoDS_Shape& aE=aItW.Value();
1453       if (aE.IsSame(theE1)) {
1454         theE2=(*(TopoDS_Edge*)(&aE));
1455         bFound=!bFound;
1456         return bFound;
1457       }
1458     }
1459   }
1460   return bFound;
1461 }
1462 //=======================================================================
1463 //function : FindFacePairs
1464 //purpose  : 
1465 //=======================================================================
1466 Standard_Boolean FindFacePairs (const TopoDS_Edge& theE,
1467                                 const BOPCol_ListOfShape& thLF,
1468                                 BOPTools_ListOfCoupleOfShape& theLCFF,
1469                                 Handle(BOPInt_Context)& theContext)
1470 {
1471   Standard_Boolean bFound;
1472   Standard_Integer i, aNbCEF;
1473   TopAbs_Orientation aOr, aOrC = TopAbs_FORWARD;
1474   BOPCol_MapOfShape aMFP;
1475   TopoDS_Face aF1, aF2;
1476   TopoDS_Edge aEL, aE1;
1477   BOPCol_ListIteratorOfListOfShape aItLF;
1478   BOPTools_CoupleOfShape aCEF, aCFF;
1479   BOPTools_ListOfCoupleOfShape aLCEF, aLCEFx;
1480   BOPTools_ListIteratorOfListOfCoupleOfShape aIt;
1481   //
1482   bFound=Standard_True;
1483   //
1484   // Preface aLCEF
1485   aItLF.Initialize(thLF);
1486   for (; aItLF.More(); aItLF.Next()) { 
1487     const TopoDS_Face& aFL=(*(TopoDS_Face*)(&aItLF.Value()));
1488     //
1489     bFound=BOPTools_AlgoTools::GetEdgeOnFace(theE, aFL, aEL);
1490     if (!bFound) {
1491       return bFound; // it can not be so
1492     }
1493     //
1494     aCEF.SetShape1(aEL);
1495     aCEF.SetShape2(aFL);
1496     aLCEF.Append(aCEF);
1497   }
1498   //
1499   aNbCEF=aLCEF.Extent();
1500   while(aNbCEF) {
1501     //
1502     // aLCEFx
1503     aLCEFx.Clear();
1504     aIt.Initialize(aLCEF);
1505     for (i=0; aIt.More(); aIt.Next(), ++i) {
1506       const BOPTools_CoupleOfShape& aCSx=aIt.Value();
1507       const TopoDS_Shape& aEx=aCSx.Shape1();
1508       const TopoDS_Shape& aFx=aCSx.Shape2();
1509       //
1510       aOr=aEx.Orientation();
1511       //
1512       if (!i) {
1513         aOrC=TopAbs::Reverse(aOr);
1514         aE1=(*(TopoDS_Edge*)(&aEx));
1515         aF1=(*(TopoDS_Face*)(&aFx));
1516         aMFP.Add(aFx);
1517         continue;
1518       }
1519       //
1520       if (aOr==aOrC) {
1521         aLCEFx.Append(aCSx);
1522         aMFP.Add(aFx);
1523       }
1524     }
1525     //
1526     // F2
1527     BOPTools_AlgoTools::GetFaceOff(aE1, aF1, aLCEFx, aF2, theContext);
1528     //
1529     aCFF.SetShape1(aF1);
1530     aCFF.SetShape2(aF2);
1531     theLCFF.Append(aCFF);
1532     //
1533     aMFP.Add(aF1);
1534     aMFP.Add(aF2);
1535     //
1536     // refine aLCEF
1537     aLCEFx.Clear();
1538     aLCEFx=aLCEF;
1539     aLCEF.Clear();
1540     aIt.Initialize(aLCEFx);
1541     for (; aIt.More(); aIt.Next()) {
1542       const BOPTools_CoupleOfShape& aCSx=aIt.Value();
1543       const TopoDS_Shape& aFx=aCSx.Shape2();
1544       if (!aMFP.Contains(aFx)) {
1545         aLCEF.Append(aCSx);
1546       }
1547     }
1548     //
1549     aNbCEF=aLCEF.Extent();
1550   }//while(aNbCEF) {
1551   //
1552   return bFound;
1553 }
1554 //=======================================================================
1555 //function : AngleWithRef
1556 //purpose  : 
1557 //=======================================================================
1558 Standard_Real AngleWithRef(const gp_Dir& theD1,
1559                            const gp_Dir& theD2,
1560                            const gp_Dir& theDRef)
1561 {
1562   Standard_Real aCosinus, aSinus, aBeta, aHalfPI, aScPr;
1563   gp_XYZ aXYZ;
1564   //
1565   aHalfPI=0.5*M_PI;
1566   //
1567   const gp_XYZ& aXYZ1=theD1.XYZ();
1568   const gp_XYZ& aXYZ2=theD2.XYZ();
1569   aXYZ=aXYZ1.Crossed(aXYZ2);
1570   aSinus=aXYZ.Modulus();
1571   aCosinus=theD1*theD2;
1572   //
1573   aBeta=0.;
1574   if (aSinus>=0.) {
1575     aBeta=aHalfPI*(1.-aCosinus);
1576   }
1577   else {
1578     aBeta=2.*M_PI-aHalfPI*(3.+aCosinus);
1579   }
1580   //
1581   aScPr=aXYZ.Dot(theDRef.XYZ());
1582   if (aScPr<0.) {
1583     aBeta=-aBeta;
1584   }
1585   return aBeta;
1586 }
1587 //=======================================================================
1588 //function : fsqrt
1589 //purpose  : 
1590 //=======================================================================
1591 Standard_Real fsqrt(Standard_Real val)  
1592 {
1593   union {
1594     int tmp;
1595     float val;
1596   } u;
1597   //
1598   u.val = (float)val;
1599   u.tmp -= 1<<23; /* Remove last bit so 1.0 gives 1.0 */
1600   /* tmp is now an approximation to logbase2(val) */
1601   u.tmp >>= 1; /* divide by 2 */
1602   u.tmp += 1<<29; /* add 64 to exponent: (e+127)/2 =(e/2)+63, */
1603   /* that represents (e/2)-64 but we want e/2 */
1604   return (double)u.val;
1605 }
1606
1607 //=======================================================================
1608 // function: IsBlockInOnFace
1609 // purpose: 
1610 //=======================================================================
1611   Standard_Boolean BOPTools_AlgoTools::IsBlockInOnFace (const IntTools_Range& aShrR,
1612                                                         const TopoDS_Face& aF,
1613                                                         const TopoDS_Edge& aE1,
1614                                                         Handle(BOPInt_Context)& aContext)
1615 {
1616   Standard_Boolean bFlag;
1617   Standard_Real f1, l1, ULD, VLD;
1618   gp_Pnt2d aP2D;
1619   gp_Pnt aP11, aP12;
1620   //
1621   aShrR.Range(f1, l1);
1622   Standard_Real dt=0.0075,  k;//dt=0.001,  k;
1623   k=dt*(l1-f1);
1624   f1=f1+k;
1625   l1=l1-k;
1626   //
1627   // Treatment P11
1628   BOPTools_AlgoTools::PointOnEdge(aE1, f1, aP11);
1629   //
1630   GeomAPI_ProjectPointOnSurf& aProjector=aContext->ProjPS(aF);
1631   aProjector.Perform(aP11);
1632   //
1633   bFlag=aProjector.IsDone();
1634   if (!bFlag) {
1635     return bFlag;
1636   }
1637   
1638   aProjector.LowerDistanceParameters(ULD, VLD);
1639   aP2D.SetCoord(ULD, VLD);
1640   //
1641   bFlag=aContext->IsPointInOnFace (aF, aP2D);
1642   //
1643   if (!bFlag) {
1644     return bFlag;
1645   }
1646   //
1647   // Treatment P12
1648   BOPTools_AlgoTools::PointOnEdge(aE1, l1, aP12);
1649   //
1650   aProjector.Perform(aP12);
1651   //
1652   bFlag=aProjector.IsDone();
1653   if (!bFlag) {
1654     return bFlag;
1655   }
1656   
1657   aProjector.LowerDistanceParameters(ULD, VLD);
1658   aP2D.SetCoord(ULD, VLD);
1659   //
1660   bFlag=aContext->IsPointInOnFace (aF, aP2D);
1661   //
1662   if (!bFlag) {
1663     return bFlag;
1664   }
1665   //
1666   // Treatment intemediate
1667   Standard_Real m1, aTolF, aTolE, aTol, aDist;
1668   m1=IntTools_Tools::IntermediatePoint(f1, l1);
1669   BOPTools_AlgoTools::PointOnEdge(aE1, m1, aP12);
1670   //
1671   aProjector.Perform(aP12);
1672   //
1673   bFlag=aProjector.IsDone();
1674   if (!bFlag) {
1675     return bFlag;
1676   }
1677   //
1678   aTolE=BRep_Tool::Tolerance(aE1);
1679   aTolF=BRep_Tool::Tolerance(aF);
1680   aTol=aTolE+aTolF;
1681   aDist=aProjector.LowerDistance();
1682   if (aDist > aTol){
1683    return Standard_False;
1684   }
1685
1686   aProjector.LowerDistanceParameters(ULD, VLD);
1687   aP2D.SetCoord(ULD, VLD);
1688   //
1689   bFlag=aContext->IsPointInOnFace (aF, aP2D);
1690   //
1691   if (!bFlag) {
1692     return bFlag;
1693   }
1694   return bFlag;
1695 }
1696
1697 //=======================================================================
1698 //function : IsMicroEdge
1699 //purpose  : 
1700 //=======================================================================
1701   Standard_Boolean BOPTools_AlgoTools::IsMicroEdge(const TopoDS_Edge& aE,
1702                                                    const Handle(BOPInt_Context)& aCtx) 
1703 {
1704   Standard_Boolean bRet;
1705   Standard_Integer iErr;
1706   Standard_Real aT1, aT2, aTmp;
1707   Handle(Geom_Curve) aC3D;
1708   TopoDS_Vertex aV1, aV2;
1709   //
1710   bRet=(BRep_Tool::Degenerated(aE) ||
1711         !BRep_Tool::IsGeometric(aE));
1712   if (bRet) {
1713     return bRet;
1714   }
1715   //
1716   aC3D=BRep_Tool::Curve(aE, aT1, aT2);
1717   TopExp::Vertices(aE, aV1, aV2);
1718   aT1=BRep_Tool::Parameter(aV1, aE);
1719   aT2=BRep_Tool::Parameter(aV2, aE);
1720   if (aT2<aT1) {
1721     aTmp=aT1;
1722     aT1=aT2;
1723     aT2=aTmp;
1724   }
1725   //
1726   BOPInt_ShrunkRange aSR;
1727   aSR.SetData(aE, aT1, aT2, aV1, aV2, aCtx);
1728   aSR.Perform();
1729   iErr=aSR.ErrorStatus();
1730   bRet = !(iErr==0);
1731   //
1732   return bRet;
1733 }
1734
1735 //=======================================================================
1736 //function : GetFaceDir
1737 //purpose  : Get binormal direction for the face in the point aP
1738 //=======================================================================
1739   void GetFaceDir(const TopoDS_Edge& aE,
1740                   const TopoDS_Face& aF,
1741                   const gp_Pnt& aP,
1742                   const Standard_Real aT,
1743                   const gp_Dir& aDTgt,
1744                   gp_Dir& aDN,
1745                   gp_Dir& aDB,
1746                   Handle(BOPInt_Context)& theContext)
1747 {
1748   BOPTools_AlgoTools3D::GetNormalToFaceOnEdge(aE, aF, aT, aDN);
1749   if (aF.Orientation()==TopAbs_REVERSED){
1750     aDN.Reverse();
1751   }
1752   aDB = aDN^aDTgt;
1753   //
1754   Handle(Geom_Surface) aS = BRep_Tool::Surface(aF);
1755   GeomAdaptor_Surface aGAS(aS);
1756   if (aGAS.GetType()!=GeomAbs_Plane) {
1757     gp_Pnt aPx;
1758     if (!FindPointInFace(aE, aF, aP, aT, aDB, aPx, theContext, aGAS)) {
1759       gp_Pnt2d aPx2D;
1760       Handle(Geom_Plane) aPL;
1761       GeomAPI_ProjectPointOnSurf aProj;
1762       //
1763       BOPTools_AlgoTools3D::PointNearEdge(aE, aF, aT, aPx2D, aPx, theContext);
1764       aPL = new Geom_Plane(aP, aDTgt);
1765       aProj.Init(aPx, aPL);
1766       aPx = aProj.NearestPoint();
1767       gp_Vec aVec(aP, aPx);
1768       aDB.SetXYZ(aVec.XYZ());
1769     }
1770   }
1771 }
1772
1773 //=======================================================================
1774 //function : FindPointInFace
1775 //purpose  : Find a point in the face in direction of <aDB>
1776 //=======================================================================
1777   Standard_Boolean FindPointInFace(const TopoDS_Edge& aE,
1778                                    const TopoDS_Face& aF,
1779                                    const gp_Pnt& aP,
1780                                    const Standard_Real /*aT*/,
1781                                    gp_Dir& aDB,
1782                                    gp_Pnt& aPOut,
1783                                    Handle(BOPInt_Context)& theContext,
1784                                    const GeomAdaptor_Surface& aGAS) 
1785 {
1786   Standard_Integer aNbItMax;
1787   Standard_Real aDt, aDtMin, aTolE, aTolF, aDist;
1788   Standard_Boolean bRet;
1789   gp_Pnt aP1;
1790   //
1791   bRet = Standard_False;
1792   aTolE = BRep_Tool::Tolerance(aE);
1793   aTolF = BRep_Tool::Tolerance(aF);
1794   aDt = 2*(aTolE+aTolF);
1795   //
1796   aDtMin=5.e-4;
1797   if (aDt < aDtMin) {
1798     Standard_Real aR;
1799     GeomAbs_SurfaceType aSType=aGAS.GetType();
1800     switch (aSType) {
1801     case GeomAbs_Cylinder:
1802       aR = aGAS.Cylinder().Radius();
1803       break;
1804     case GeomAbs_Cone:
1805       aR = aGAS.Cone().RefRadius();
1806       break;
1807     case GeomAbs_Sphere:
1808       aR = aGAS.Sphere().Radius();
1809       break;
1810     case GeomAbs_Torus:
1811       aR = aGAS.Torus().MinorRadius();
1812       break;
1813     case GeomAbs_SurfaceOfRevolution:
1814       aR=1.;
1815       break;
1816     default:
1817       aR=0.001;
1818     }
1819     if (aR < 0.01) {
1820       aDtMin=5.e-6;
1821     }
1822     else if (aR > 100.) {
1823       aDtMin*=10;
1824     }
1825     if (aDt < aDtMin) {
1826       aDt=aDtMin;
1827     }
1828   }
1829   //
1830   GeomAPI_ProjectPointOnSurf& aProj=theContext->ProjPS(aF);
1831   aNbItMax = 15;
1832   //
1833   do {
1834     aP1.SetCoord(aP.X()+aDt*aDB.X(),
1835                  aP.Y()+aDt*aDB.Y(),
1836                  aP.Z()+aDt*aDB.Z());
1837     //
1838     aProj.Perform(aP1);
1839     if (!aProj.IsDone()) {
1840       return bRet;
1841     }
1842     aPOut = aProj.NearestPoint();
1843     aDist = aProj.LowerDistance();
1844     //
1845     gp_Vec aV(aP, aPOut);
1846     aDB.SetXYZ(aV.XYZ());
1847   } while (aDist>Precision::Angular() && --aNbItMax);
1848   //
1849   bRet = aDist < Precision::Angular();
1850   return bRet;
1851 }