0032539: Modeling Algorithms - Parallelize BRepExtrema_DistShapeShape algorithm
[occt.git] / src / BRepExtrema / BRepExtrema_DistShapeShape.cxx
1 // Created on: 1996-04-22
2 // Created by: Herve LOUESSARD
3 // Copyright (c) 1996-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and/or modify it under
9 // the terms of the GNU Lesser General Public License version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 // Modified: Mps(10-04-97) portage WNT
18
19 #include <BRepExtrema_DistShapeShape.hxx>
20
21 #include <Standard_OStream.hxx>
22 #include <TopTools_IndexedMapOfShape.hxx>
23 #include <BRepBndLib.hxx>
24 #include <Bnd_Box.hxx>
25 #include <TopExp.hxx>
26 #include <BRepExtrema_DistanceSS.hxx>
27 #include <TopoDS.hxx>
28 #include <TopAbs.hxx>
29 #include <TopoDS_Vertex.hxx>
30 #include <TopoDS_Edge.hxx>
31 #include <TopoDS_Face.hxx>
32 #include <TopAbs_ShapeEnum.hxx>
33 #include <Precision.hxx>
34 #include <Bnd_SeqOfBox.hxx>
35 #include <BRepExtrema_UnCompatibleShape.hxx>
36 #include <BRep_Tool.hxx>
37 #include <BRepClass3d_SolidClassifier.hxx>
38 #include <NCollection_Vector.hxx>
39 #include <OSD_Parallel.hxx>
40 #include <StdFail_NotDone.hxx>
41
42 #include <algorithm>
43 namespace
44 {
45
46   static void Decomposition(const TopoDS_Shape&         S,
47                             TopTools_IndexedMapOfShape& MapV,
48                             TopTools_IndexedMapOfShape& MapE,
49                             TopTools_IndexedMapOfShape& MapF)
50   {
51     MapV.Clear();
52     MapE.Clear();
53     MapF.Clear();
54     TopExp::MapShapes(S,TopAbs_VERTEX,MapV);
55     TopExp::MapShapes(S,TopAbs_EDGE,MapE);
56     TopExp::MapShapes(S,TopAbs_FACE,MapF);
57   }
58
59   static void BoxCalculation(const TopTools_IndexedMapOfShape& Map,
60                              Bnd_Array1OfBox&                  SBox)
61   {
62     for (Standard_Integer i = 1; i <= Map.Extent(); i++)
63     {
64       Bnd_Box box;
65       BRepBndLib::Add(Map(i), box);
66       SBox[i] = box;
67     }
68   }
69
70   inline Standard_Real DistanceInitiale(const TopoDS_Vertex V1,
71                                         const TopoDS_Vertex V2)
72   {
73     return (BRep_Tool::Pnt(V1).Distance(BRep_Tool::Pnt(V2)));
74   }
75
76   //! Pair of objects to check extrema.
77   struct BRepExtrema_CheckPair
78   {
79     Standard_Integer Index1;   //!< Index of the 1st sub-shape
80     Standard_Integer Index2;   //!< Index of the 2nd sub-shape
81     Standard_Real    Distance; //!< Distance between sub-shapes
82
83     //! Uninitialized constructor for collection.
84     BRepExtrema_CheckPair()
85     : Index1(0),
86       Index2(0),
87       Distance(0.0)
88     {
89     }
90
91     //! Creates new pair of sub-shapes.
92     BRepExtrema_CheckPair (Standard_Integer theIndex1,
93                            Standard_Integer theIndex2,
94                            Standard_Real    theDistance)
95     : Index1   (theIndex1),
96       Index2   (theIndex2),
97       Distance (theDistance) {}
98   };
99
100   // Used by std::sort function
101   static Standard_Boolean BRepExtrema_CheckPair_Comparator (const BRepExtrema_CheckPair& theLeft,
102                                                             const BRepExtrema_CheckPair& theRight)
103   {
104     return (theLeft.Distance < theRight.Distance);
105   }
106 }
107
108 //=======================================================================
109 //struct   : IndexBand
110 //purpose  : 
111 //=======================================================================
112 struct IndexBand
113 {
114   IndexBand():
115     First(0),
116     Last(0)
117   {
118   }  
119   
120   IndexBand(Standard_Integer theFirtsIndex,
121             Standard_Integer theLastIndex):
122     First(theFirtsIndex),
123     Last(theLastIndex)
124   {
125   }
126   Standard_Integer First;
127   Standard_Integer Last;
128 };
129
130 //=======================================================================
131 //struct   : ThreadSolution
132 //purpose  : 
133 //=======================================================================
134 struct ThreadSolution
135 {
136   ThreadSolution(Standard_Integer theTaskNum):
137     Shape1(0, theTaskNum-1),
138     Shape2(0, theTaskNum-1),
139     Dist(0, theTaskNum-1)
140   {
141     Dist.Init(DBL_MAX);
142   }
143
144   NCollection_Array1<BRepExtrema_SeqOfSolution> Shape1;
145   NCollection_Array1<BRepExtrema_SeqOfSolution> Shape2;
146   NCollection_Array1<Standard_Real>             Dist;
147 };
148
149 //=======================================================================
150 //struct   : VertexFunctor
151 //purpose  : 
152 //=======================================================================
153 struct VertexFunctor
154 {
155   VertexFunctor(NCollection_Array1<IndexBand>* theBandArray,
156                 const Message_ProgressRange& theRange):
157     BandArray(theBandArray),
158     Solution(theBandArray->Size()),
159     Map1(NULL),
160     Map2(NULL),
161     Scope(theRange, "Vertices distances calculating", theBandArray->Size()),
162     Ranges(0, theBandArray->Size() - 1),
163     Eps(Precision::Confusion()),
164     StartDist(0.0)
165   {
166     for (Standard_Integer i = 0; i < theBandArray->Size(); ++i)
167     {
168       Ranges.SetValue(i, Scope.Next());
169     }
170   }
171
172   void operator() (const Standard_Integer theIndex) const
173   {
174     const Standard_Integer aCount2 = Map2->Extent();
175     const Standard_Integer aFirst = BandArray->Value(theIndex).First;
176     const Standard_Integer aLast  = BandArray->Value(theIndex).Last;
177     Solution.Dist[theIndex] = StartDist;
178
179     Message_ProgressScope aScope(Ranges[theIndex], NULL, (double)aLast - aFirst);
180
181
182     for (Standard_Integer anIdx1 = aFirst; anIdx1 <= aLast; ++anIdx1)
183     {
184       if (!aScope.More())
185       {
186         break;
187       }
188       aScope.Next();
189
190       const TopoDS_Vertex& aVertex1 = TopoDS::Vertex(Map1->FindKey(anIdx1));
191       const gp_Pnt aPoint1 = BRep_Tool::Pnt(aVertex1);
192       for (Standard_Integer anIdx2 = 1; anIdx2 <= aCount2; ++anIdx2)
193       {
194         const TopoDS_Vertex& aVertex2 = TopoDS::Vertex(Map2->FindKey(anIdx2));
195         const gp_Pnt aPoint2 = BRep_Tool::Pnt(aVertex2);
196         const Standard_Real aDist = aPoint1.Distance(aPoint2);
197         {
198           if (aDist < Solution.Dist[theIndex] - Eps)
199           { 
200             const BRepExtrema_SolutionElem Sol1(aDist, aPoint1, BRepExtrema_IsVertex, aVertex1);
201             const BRepExtrema_SolutionElem Sol2(aDist, aPoint2, BRepExtrema_IsVertex, aVertex2);
202
203             Solution.Shape1[theIndex].Clear();
204             Solution.Shape2[theIndex].Clear();
205             Solution.Shape1[theIndex].Append(Sol1);
206             Solution.Shape2[theIndex].Append(Sol2);
207             
208             Solution.Dist[theIndex] = aDist;
209           }
210           else if (Abs(aDist - Solution.Dist[theIndex]) < Eps)
211           {
212             const BRepExtrema_SolutionElem Sol1(aDist, aPoint1, BRepExtrema_IsVertex, aVertex1);
213             const BRepExtrema_SolutionElem Sol2(aDist, aPoint2, BRepExtrema_IsVertex, aVertex2);
214             Solution.Shape1[theIndex].Append(Sol1);
215             Solution.Shape2[theIndex].Append(Sol2);
216
217             if (Solution.Dist[theIndex] > aDist)
218             {
219               Solution.Dist[theIndex] = aDist;
220             }
221           }          
222         }
223       }
224     }
225   }
226
227   NCollection_Array1<IndexBand>*            BandArray;
228   mutable ThreadSolution                    Solution;
229   const TopTools_IndexedMapOfShape*         Map1;
230   const TopTools_IndexedMapOfShape*         Map2;
231   Message_ProgressScope                     Scope;
232   NCollection_Array1<Message_ProgressRange> Ranges;
233   Standard_Real                             Eps;
234   Standard_Real                             StartDist;
235 };
236
237 //=======================================================================
238 //function : DistanceVertVert
239 //purpose  : 
240 //=======================================================================
241 Standard_Boolean BRepExtrema_DistShapeShape::DistanceVertVert(const TopTools_IndexedMapOfShape& theMap1,
242                                                               const TopTools_IndexedMapOfShape& theMap2,
243                                                               const Message_ProgressRange& theRange)
244 {
245   const Standard_Integer aCount1 = theMap1.Extent();
246   const Standard_Integer aMinTaskSize = aCount1 < 10 ? aCount1 : 10;
247   const Handle(OSD_ThreadPool)& aThreadPool = OSD_ThreadPool::DefaultPool();
248   const Standard_Integer aNbThreads = aThreadPool->NbThreads();
249   Standard_Integer aNbTasks = aNbThreads;
250   Standard_Integer aTaskSize = (Standard_Integer) Ceiling((double) aCount1 / aNbTasks);
251   if (aTaskSize < aMinTaskSize)
252   {
253     aTaskSize = aMinTaskSize;
254     aNbTasks = (Standard_Integer) Ceiling((double) aCount1 / aTaskSize);
255   }
256
257   Standard_Integer aFirstIndex(1);
258   NCollection_Array1<IndexBand> aBandArray(0, aNbTasks - 1);
259   Message_ProgressScope aDistScope(theRange, NULL, 1);
260
261   for (Standard_Integer anI = 0; anI < aBandArray.Size(); ++anI)
262   {
263     if (aCount1 < aFirstIndex + aTaskSize - 1)
264     {
265       aTaskSize = aCount1 - aFirstIndex + 1;
266     }
267     aBandArray.SetValue(anI, IndexBand(aFirstIndex, aFirstIndex + aTaskSize - 1));
268     aFirstIndex += aTaskSize;
269   }
270
271   VertexFunctor aFunctor(&aBandArray, aDistScope.Next());
272   aFunctor.Map1            = &theMap1;
273   aFunctor.Map2            = &theMap2;
274   aFunctor.StartDist       = myDistRef;
275   aFunctor.Eps             = myEps;
276
277   OSD_Parallel::For(0, aNbTasks, aFunctor, !myIsMultiThread);
278   if (!aDistScope.More())
279   {
280     return Standard_False;
281   }    
282   for (Standard_Integer anI = 0; anI < aFunctor.Solution.Dist.Size(); ++anI)
283   {
284     Standard_Real aDist = aFunctor.Solution.Dist[anI];
285     if (aDist < myDistRef - myEps)
286     {
287       mySolutionsShape1.Clear();
288       mySolutionsShape2.Clear();
289       mySolutionsShape1.Append(aFunctor.Solution.Shape1[anI]);
290       mySolutionsShape2.Append(aFunctor.Solution.Shape2[anI]);
291       myDistRef = aDist;
292     }
293     else if (Abs(aDist - myDistRef) < myEps)
294     {
295       mySolutionsShape1.Append(aFunctor.Solution.Shape1[anI]);
296       mySolutionsShape2.Append(aFunctor.Solution.Shape2[anI]);
297       myDistRef = aDist;
298     }
299   }
300   return Standard_True;
301 }
302
303 //=======================================================================
304 //struct   : DistanceFunctor
305 //purpose  : 
306 //=======================================================================
307 struct DistanceFunctor
308 {
309   DistanceFunctor(NCollection_Array1<NCollection_Array1<BRepExtrema_CheckPair> >* theArrayOfArrays,
310                   const Message_ProgressRange& theRange):
311     ArrayOfArrays(theArrayOfArrays), 
312     Solution(ArrayOfArrays->Size()),
313     Map1(NULL),
314     Map2(NULL),
315     LBox1(NULL),
316     LBox2(NULL),
317     Scope(theRange, "Shapes distances calculating", theArrayOfArrays->Size()),
318     Ranges(0, theArrayOfArrays->Size() - 1),
319     Eps(Precision::Confusion()),
320     StartDist(0.0)
321   {
322     for (Standard_Integer i = 0; i < theArrayOfArrays->Size(); ++i)
323     {
324       Ranges.SetValue(i, Scope.Next());
325     }
326   }
327
328   void operator() (const Standard_Integer theIndex) const
329   {
330     Message_ProgressScope aScope(Ranges[theIndex], NULL, ArrayOfArrays->Value(theIndex).Size());
331     Solution.Dist[theIndex] = StartDist;
332     for (Standard_Integer i = 0; i < ArrayOfArrays->Value(theIndex).Size(); i++)
333     {
334       if (!aScope.More())
335       {
336         return;
337       }
338       aScope.Next();
339       const BRepExtrema_CheckPair& aPair = ArrayOfArrays->Value(theIndex).Value(i);
340       if (aPair.Distance > Solution.Dist[theIndex] + Eps)
341       {
342         break; // early search termination
343       }
344       const Bnd_Box& aBox1 = LBox1->Value(aPair.Index1);
345       const Bnd_Box& aBox2 = LBox2->Value(aPair.Index2);
346       const TopoDS_Shape& aShape1 = Map1->FindKey(aPair.Index1);
347       const TopoDS_Shape& aShape2 = Map2->FindKey(aPair.Index2);
348       BRepExtrema_DistanceSS aDistTool(aShape1, aShape2, aBox1, aBox2, Solution.Dist[theIndex], Eps);
349       const Standard_Real aDist = aDistTool.DistValue();
350       if (aDistTool.IsDone())
351       {
352         if (aDist < Solution.Dist[theIndex] - Eps)
353         {
354           Solution.Shape1[theIndex].Clear();
355           Solution.Shape2[theIndex].Clear();
356
357           BRepExtrema_SeqOfSolution aSeq1 = aDistTool.Seq1Value();
358           BRepExtrema_SeqOfSolution aSeq2 = aDistTool.Seq2Value();
359
360           Solution.Shape1[theIndex].Append(aSeq1);
361           Solution.Shape2[theIndex].Append(aSeq2);
362
363           Solution.Dist[theIndex] = aDistTool.DistValue();
364         }
365         else if (Abs(aDist - Solution.Dist[theIndex]) < Eps)
366         {
367           BRepExtrema_SeqOfSolution aSeq1 = aDistTool.Seq1Value();
368           BRepExtrema_SeqOfSolution aSeq2 = aDistTool.Seq2Value();
369
370           Solution.Shape1[theIndex].Append(aSeq1);
371           Solution.Shape2[theIndex].Append(aSeq2);
372           if (Solution.Dist[theIndex] > aDist)
373           {
374             Solution.Dist[theIndex] = aDist;
375           }
376         }
377       }
378     }
379   }
380
381   NCollection_Array1<NCollection_Array1<BRepExtrema_CheckPair> >* ArrayOfArrays;
382   mutable ThreadSolution                    Solution;
383   const TopTools_IndexedMapOfShape*         Map1;
384   const TopTools_IndexedMapOfShape*         Map2;
385   const Bnd_Array1OfBox*                    LBox1;
386   const Bnd_Array1OfBox*                    LBox2;
387   Message_ProgressScope                     Scope;
388   NCollection_Array1<Message_ProgressRange> Ranges;
389   Standard_Real                             Eps;
390   Standard_Real                             StartDist;
391 };
392
393
394 //=======================================================================
395 //struct   : DistancePairFunctor
396 //purpose  : 
397 //=======================================================================
398 struct DistancePairFunctor
399 {
400   DistancePairFunctor(NCollection_Array1<IndexBand>* theBandArray,
401                       const Message_ProgressRange& theRange):
402     BandArray(theBandArray),
403     PairList(0, theBandArray->Size() - 1),
404     LBox1(NULL),
405     LBox2(NULL),
406     Scope(theRange, "Boxes distances calculating", theBandArray->Size()),
407     Ranges(0, theBandArray->Size() - 1),
408     DistRef(0),
409     Eps(Precision::Confusion())
410   {
411     for (Standard_Integer i = 0; i < theBandArray->Size(); ++i)
412     {
413       Ranges.SetValue(i, Scope.Next());
414     }
415   }
416
417   void operator() (const Standard_Integer theIndex) const
418   {
419     const Standard_Integer aFirst = BandArray->Value(theIndex).First;
420     const Standard_Integer aLast  = BandArray->Value(theIndex).Last;
421
422     Message_ProgressScope aScope(Ranges[theIndex], NULL, (double) aLast - aFirst);
423
424     for (Standard_Integer anIdx1 = aFirst; anIdx1 <= aLast; ++anIdx1)
425     {
426       if (!aScope.More())
427       {
428         break;
429       }
430       aScope.Next();
431
432       for (Standard_Integer anIdx2 = 1; anIdx2 <= LBox2->Size(); ++anIdx2)
433       {
434         const Bnd_Box& aBox1 = LBox1->Value(anIdx1);
435         const Bnd_Box& aBox2 = LBox2->Value(anIdx2);
436         if (aBox1.IsVoid() || aBox2.IsVoid())
437         {
438           continue;
439         }
440
441         const Standard_Real aDist = aBox1.Distance(aBox2);       
442         if (aDist - DistRef < Eps)
443         {
444           PairList[theIndex].Append(BRepExtrema_CheckPair(anIdx1, anIdx2, aDist));
445         }
446       }
447     }
448   }
449
450   Standard_Integer ListSize()
451   {
452     Standard_Integer aSize(0);
453     for (Standard_Integer anI = PairList.Lower(); anI <= PairList.Upper(); ++anI)
454     {
455       aSize += PairList[anI].Size();
456     }
457     return aSize;
458   }
459
460   NCollection_Array1<IndexBand>*             BandArray;
461   mutable NCollection_Array1<NCollection_Vector<BRepExtrema_CheckPair> > PairList;
462   const Bnd_Array1OfBox*                     LBox1;
463   const Bnd_Array1OfBox*                     LBox2;
464   Message_ProgressScope                      Scope;
465   NCollection_Array1<Message_ProgressRange>  Ranges;
466   Standard_Real                              DistRef;
467   Standard_Real                              Eps;
468 };
469
470 //=======================================================================
471 //function : DistanceMapMap
472 //purpose  : 
473 //=======================================================================
474
475 Standard_Boolean BRepExtrema_DistShapeShape::DistanceMapMap (const TopTools_IndexedMapOfShape& theMap1,
476                                                              const TopTools_IndexedMapOfShape& theMap2,
477                                                              const Bnd_Array1OfBox&            theLBox1,
478                                                              const Bnd_Array1OfBox&            theLBox2,
479                                                              const Message_ProgressRange&      theRange)
480 {
481   const Standard_Integer aCount1 = theMap1.Extent();
482   const Standard_Integer aCount2 = theMap2.Extent();
483
484   if (aCount1 == 0 || aCount2 == 0)
485   {
486     return Standard_True;
487   }
488
489   Message_ProgressScope aTwinScope(theRange, NULL, 1.0);
490
491   const Handle(OSD_ThreadPool)& aThreadPool = OSD_ThreadPool::DefaultPool();
492   const Standard_Integer aNbThreads = aThreadPool->NbThreads();
493   const Standard_Integer aMinPairTaskSize = aCount1 < 10 ? aCount1 : 10;
494   Standard_Integer aNbPairTasks = aNbThreads;
495   Standard_Integer aPairTaskSize = (Standard_Integer) Ceiling((double) aCount1 / aNbPairTasks);
496   if (aPairTaskSize < aMinPairTaskSize)
497   {
498     aPairTaskSize = aMinPairTaskSize;
499     aNbPairTasks = (Standard_Integer) Ceiling((double) aCount1 / aPairTaskSize);
500   }
501
502   Standard_Integer aFirstIndex(1);
503   NCollection_Array1<IndexBand> aBandArray(0, aNbPairTasks - 1);
504
505   for (Standard_Integer anI = 0; anI < aBandArray.Size(); ++anI)
506   {
507     if (aCount1 < aFirstIndex + aPairTaskSize - 1)
508     {
509       aPairTaskSize = aCount1 - aFirstIndex + 1;
510     }
511     aBandArray.SetValue(anI, IndexBand(aFirstIndex, aFirstIndex + aPairTaskSize - 1));
512     aFirstIndex += aPairTaskSize;
513   }
514
515   aTwinScope.Next(0.15);
516   DistancePairFunctor aPairFunctor(&aBandArray, aTwinScope.Next(0.15));
517   aPairFunctor.LBox1 = &theLBox1;
518   aPairFunctor.LBox2 = &theLBox2;
519   aPairFunctor.DistRef = myDistRef;
520   aPairFunctor.Eps = myEps;
521
522   OSD_Parallel::For(0, aNbPairTasks, aPairFunctor, !myIsMultiThread);
523   if (!aTwinScope.More())
524   {
525     return Standard_False;
526   }
527   Standard_Integer aListSize = aPairFunctor.ListSize();
528   if(aListSize == 0)
529   {
530     return Standard_True;
531   }
532   NCollection_Array1<BRepExtrema_CheckPair> aPairList(0, aListSize-1);
533   Standard_Integer aListIndex(0);
534   for (Standard_Integer anI = 0; anI < aPairFunctor.PairList.Size(); ++anI)
535   {
536     for (Standard_Integer aJ = 0; aJ < aPairFunctor.PairList[anI].Size(); ++aJ)
537     {
538       aPairList[aListIndex] = aPairFunctor.PairList[anI][aJ];
539       ++aListIndex;
540     }
541   }
542
543   std::stable_sort(aPairList.begin(), aPairList.end(), BRepExtrema_CheckPair_Comparator);
544
545   const Standard_Integer aMapSize = aPairList.Size();
546   Standard_Integer aNbTasks = aMapSize < aNbThreads ? aMapSize : aNbThreads;
547   Standard_Integer aTaskSize = (Standard_Integer) Ceiling((double) aMapSize / aNbTasks);
548  
549   NCollection_Array1<NCollection_Array1<BRepExtrema_CheckPair> > anArrayOfArray(0, aNbTasks - 1);
550   // Since aPairList is sorted in ascending order of distances between Bnd_Boxes,
551   // BRepExtrema_CheckPair are distributed to tasks one by one from smallest to largest,
552   // and not ranges, as for DistancePairFunctor.
553   // Since aMapSize may not be divisible entirely by the number of tasks,
554   // some tasks should receive one BRepExtrema_CheckPair less than the rest.
555   // aLastRowLimit defines the task number from which to start tasks containing
556   // fewer BRepExtrema_CheckPair
557   Standard_Integer aLastRowLimit = ((aMapSize % aNbTasks) == 0) ? aNbTasks : (aMapSize % aNbTasks);
558   for (Standard_Integer anI = 0; anI < aTaskSize; ++anI)
559   {
560     for (Standard_Integer aJ = 0; aJ < aNbTasks; ++aJ)
561     {
562       if (anI == 0)
563       {
564         Standard_Integer aVectorSize = aTaskSize;
565         if (aJ >= aLastRowLimit)
566         {
567           aVectorSize--;
568         }
569         anArrayOfArray[aJ].Resize(0, aVectorSize - 1, Standard_False);
570       }
571       if (anI < anArrayOfArray[aJ].Size())
572       {
573         anArrayOfArray[aJ][anI] = aPairList(anI*aNbTasks + aJ);
574       }
575       else
576       {
577         break;
578       }
579     }
580   }
581   DistanceFunctor aFunctor(&anArrayOfArray, aTwinScope.Next(0.85));
582   aFunctor.Map1 = &theMap1;
583   aFunctor.Map2 = &theMap2;
584   aFunctor.LBox1 = &theLBox1;
585   aFunctor.LBox2 = &theLBox2;
586   aFunctor.Eps = myEps;
587   aFunctor.StartDist = myDistRef;
588
589   OSD_Parallel::For(0, aNbTasks, aFunctor, !myIsMultiThread);
590   if (!aTwinScope.More())
591   {
592     return Standard_False;
593   }
594
595   for (Standard_Integer anI = 0; anI < aFunctor.Solution.Dist.Size(); ++anI)
596   {
597     Standard_Real aDist = aFunctor.Solution.Dist[anI];
598     if (aDist < myDistRef - myEps)
599     {
600       mySolutionsShape1.Clear();
601       mySolutionsShape2.Clear();
602       mySolutionsShape1.Append(aFunctor.Solution.Shape1[anI]);
603       mySolutionsShape2.Append(aFunctor.Solution.Shape2[anI]);
604       myDistRef = aDist;
605     }
606     else if (Abs(aDist - myDistRef) < myEps)
607     {
608       mySolutionsShape1.Append(aFunctor.Solution.Shape1[anI]);
609       mySolutionsShape2.Append(aFunctor.Solution.Shape2[anI]);
610       if (myDistRef > aDist)
611       {
612         myDistRef = aDist;
613       }
614     }
615   }
616   return Standard_True;
617 }
618
619 //=======================================================================
620 //function : BRepExtrema_DistShapeShape
621 //purpose  : 
622 //=======================================================================
623
624 BRepExtrema_DistShapeShape::BRepExtrema_DistShapeShape()
625 : myDistRef (0.0),
626   myIsDone (Standard_False),
627   myInnerSol (Standard_False),
628   myEps (Precision::Confusion()),
629   myIsInitS1 (Standard_False),
630   myIsInitS2 (Standard_False),
631   myFlag (Extrema_ExtFlag_MINMAX),
632   myAlgo (Extrema_ExtAlgo_Grad),
633   myIsMultiThread(Standard_False)
634 {
635 }
636
637 //=======================================================================
638 //function : BRepExtrema_DistShapeShape
639 //purpose  : 
640 //=======================================================================
641 BRepExtrema_DistShapeShape::BRepExtrema_DistShapeShape(const TopoDS_Shape& Shape1,
642                                                        const TopoDS_Shape& Shape2,
643                                                        const Extrema_ExtFlag F,
644                                                        const Extrema_ExtAlgo A,
645                                                        const Message_ProgressRange& theRange)
646 : myDistRef (0.0),
647   myIsDone (Standard_False),
648   myInnerSol (Standard_False),
649   myEps (Precision::Confusion()),
650   myIsInitS1 (Standard_False),
651   myIsInitS2 (Standard_False),
652   myFlag (F),
653   myAlgo (A),
654   myIsMultiThread(Standard_False)
655 {
656   LoadS1(Shape1);
657   LoadS2(Shape2);
658   Perform(theRange);
659 }
660
661 //=======================================================================
662 //function : BRepExtrema_DistShapeShape
663 //purpose  : 
664 //=======================================================================
665
666 BRepExtrema_DistShapeShape::BRepExtrema_DistShapeShape(const TopoDS_Shape& Shape1,
667                                                        const TopoDS_Shape& Shape2,
668                                                        const Standard_Real theDeflection,
669                                                        const Extrema_ExtFlag F,
670                                                        const Extrema_ExtAlgo A,
671                                                        const Message_ProgressRange& theRange)
672 : myDistRef (0.0),
673   myIsDone (Standard_False),
674   myInnerSol (Standard_False),
675   myEps (theDeflection),
676   myIsInitS1 (Standard_False),
677   myIsInitS2 (Standard_False),
678   myFlag (F),
679   myAlgo (A),
680   myIsMultiThread(Standard_False)
681 {
682   LoadS1(Shape1);
683   LoadS2(Shape2);
684   Perform(theRange);
685 }
686
687 //=======================================================================
688 //function : LoadS1
689 //purpose  : 
690 //=======================================================================
691
692 void BRepExtrema_DistShapeShape::LoadS1 (const TopoDS_Shape& Shape1)
693 {
694   myShape1 = Shape1;
695   myIsInitS1 = Standard_False;
696   Decomposition (Shape1, myMapV1, myMapE1, myMapF1);
697 }
698
699 //=======================================================================
700 //function : LoadS2
701 //purpose  : 
702 //=======================================================================
703
704 void BRepExtrema_DistShapeShape::LoadS2 (const TopoDS_Shape& Shape2)
705 {
706   myShape2 = Shape2;
707   myIsInitS2 = Standard_False;
708   Decomposition (Shape2, myMapV2, myMapE2, myMapF2);
709 }
710
711 //=======================================================================
712 //struct   : TreatmentFunctor
713 //purpose  : 
714 //=======================================================================
715 struct TreatmentFunctor
716 {
717   TreatmentFunctor(NCollection_Array1<NCollection_Array1<TopoDS_Shape> >* theArrayOfArrays,
718                    const Message_ProgressRange& theRange):
719   ArrayOfArrays(theArrayOfArrays),
720   SolutionsShape1(NULL),
721   SolutionsShape2(NULL),
722   Scope(theRange, "Search for the inner solid", theArrayOfArrays->Size()),
723   Ranges(0, theArrayOfArrays->Size() - 1),
724   DistRef(0),
725   InnerSol(NULL),
726   IsDone(NULL),
727   Mutex(NULL)
728   {
729     for (Standard_Integer i = 0; i < theArrayOfArrays->Size(); ++i)
730     {
731       Ranges.SetValue(i, Scope.Next());
732     }
733   }
734
735   void operator() (const Standard_Integer theIndex) const
736   {
737     const Standard_Real aTolerance = 0.001;
738     Message_ProgressScope aScope(Ranges[theIndex], NULL, ArrayOfArrays->Value(theIndex).Size());
739     BRepClass3d_SolidClassifier aClassifier(Shape);
740
741     for (Standard_Integer i = 0; i < ArrayOfArrays->Value(theIndex).Size(); i++)
742     {
743       if (!aScope.More())
744       {
745         break;
746       }
747       aScope.Next();
748       if (*IsDone)
749       {
750         break;
751       }
752
753       const TopoDS_Vertex& aVertex = TopoDS::Vertex(ArrayOfArrays->Value(theIndex).Value(i));
754       const gp_Pnt aPnt = BRep_Tool::Pnt(aVertex);
755       aClassifier.Perform(aPnt, aTolerance);
756       if (aClassifier.State() == TopAbs_IN)
757       {
758         Standard_Mutex::Sentry aLock(Mutex.get());
759         *InnerSol = Standard_True;
760         *DistRef  = 0.;
761         *IsDone   = Standard_True;
762         BRepExtrema_SolutionElem aSolElem(0, aPnt, BRepExtrema_IsVertex, aVertex);
763         SolutionsShape1->Append(aSolElem);
764         SolutionsShape2->Append(aSolElem);
765         break;
766       }
767     }
768   }
769
770   NCollection_Array1<NCollection_Array1<TopoDS_Shape> >* ArrayOfArrays;
771   BRepExtrema_SeqOfSolution*                 SolutionsShape1;
772   BRepExtrema_SeqOfSolution*                 SolutionsShape2;
773   TopoDS_Shape                               Shape;
774   Message_ProgressScope                      Scope;
775   NCollection_Array1<Message_ProgressRange>  Ranges;
776   Standard_Real*                             DistRef;
777   volatile Standard_Boolean*                 InnerSol;
778   volatile Standard_Boolean*                 IsDone;
779   Handle(Standard_HMutex)                    Mutex;
780 };
781
782 //=======================================================================
783 //function : SolidTreatment
784 //purpose  : 
785 //=======================================================================
786 Standard_Boolean BRepExtrema_DistShapeShape::SolidTreatment(const TopoDS_Shape& theShape,
787                                                             const TopTools_IndexedMapOfShape& theVertexMap,
788                                                             const Message_ProgressRange& theRange)
789 {
790   const Standard_Integer aMapSize = theVertexMap.Extent();
791   const Standard_Integer aMinTaskSize = 3;
792   const Handle(OSD_ThreadPool)& aThreadPool = OSD_ThreadPool::DefaultPool();
793   const Standard_Integer aNbThreads = aThreadPool->NbThreads();
794   Standard_Integer aNbTasks = aNbThreads * 10;
795   Standard_Integer aTaskSize = (Standard_Integer) Ceiling((double) aMapSize / aNbTasks);
796   if (aTaskSize < aMinTaskSize)
797   {
798     aTaskSize = aMinTaskSize;
799     aNbTasks = (Standard_Integer) Ceiling((double) aMapSize / aTaskSize);
800   }
801
802   NCollection_Array1< NCollection_Array1<TopoDS_Shape> > anArrayOfArray(0, aNbTasks - 1);
803   for (Standard_Integer anI = 1; anI <= aMapSize; ++anI)
804   {
805     Standard_Integer aVectIndex = (anI - 1) / aTaskSize;
806     Standard_Integer aShapeIndex = (anI - 1) % aTaskSize;
807     if (aShapeIndex == 0)
808     {
809       Standard_Integer aVectorSize = aTaskSize;
810       Standard_Integer aTailSize = aMapSize - aVectIndex * aTaskSize;
811       if (aTailSize < aTaskSize)
812       {
813         aVectorSize = aTailSize;
814       }
815       anArrayOfArray[aVectIndex].Resize(0, aVectorSize - 1, Standard_False);
816     }
817     anArrayOfArray[aVectIndex][aShapeIndex] = theVertexMap(anI);
818   }
819
820   Message_ProgressScope aScope(theRange, "Solid treatment", aNbTasks);
821   TreatmentFunctor aFunctor(&anArrayOfArray, aScope.Next());
822   aFunctor.SolutionsShape1 = &mySolutionsShape1;
823   aFunctor.SolutionsShape2 = &mySolutionsShape2;
824   aFunctor.Shape = theShape;
825   aFunctor.DistRef = &myDistRef;
826   aFunctor.InnerSol = &myInnerSol;
827   aFunctor.IsDone = &myIsDone;
828   if (myIsMultiThread)
829   {
830     aFunctor.Mutex.reset(new Standard_HMutex());
831   }
832
833   OSD_Parallel::For(0, aNbTasks, aFunctor, !myIsMultiThread);
834
835   if (!aScope.More())
836   {
837     return Standard_False;
838   }
839   return Standard_True;
840 }
841
842 //=======================================================================
843 //function : Perform
844 //purpose  : 
845 //=======================================================================
846
847 Standard_Boolean BRepExtrema_DistShapeShape::Perform(const Message_ProgressRange& theRange)
848 {
849   myIsDone = Standard_False;
850   myInnerSol = Standard_False;
851   mySolutionsShape1.Clear();
852   mySolutionsShape2.Clear();
853
854   if (myShape1.IsNull() || myShape2.IsNull())
855     return Standard_False;
856
857   // Treatment of solids
858   Standard_Boolean anIsSolid1 = (myShape1.ShapeType() == TopAbs_SOLID) ||
859     (myShape1.ShapeType() == TopAbs_COMPSOLID);
860   Standard_Boolean anIsSolid2 = (myShape2.ShapeType() == TopAbs_SOLID) ||
861     (myShape2.ShapeType() == TopAbs_COMPSOLID);
862   Standard_Integer aRootStepsNum = 9; // By num of DistanceMapMap calls
863   aRootStepsNum = anIsSolid1 ? aRootStepsNum + 1 : aRootStepsNum;
864   aRootStepsNum = anIsSolid2 ? aRootStepsNum + 1 : aRootStepsNum;
865   Message_ProgressScope aRootScope(theRange, "calculating distance", aRootStepsNum);
866
867   if (anIsSolid1)
868   {
869     if (!SolidTreatment(myShape1, myMapV2, aRootScope.Next()))
870     {
871       return Standard_False;
872     }
873   }
874
875   if (anIsSolid2 && (!myInnerSol))
876   {
877     if (!SolidTreatment(myShape2, myMapV1, aRootScope.Next()))
878     {
879       return Standard_False;
880     }
881   }
882
883   if (!myInnerSol)
884   {
885     if (!myIsInitS1) // rebuild cached data for 1st shape
886     {
887       if (!myMapV1.IsEmpty())
888       {
889         myBV1.Resize(1, myMapV1.Extent(), Standard_False);
890       }
891       if (!myMapE1.IsEmpty())
892       {
893         myBE1.Resize(1, myMapE1.Extent(), Standard_False);
894       }
895       if (!myMapF1.IsEmpty())
896       {
897         myBF1.Resize(1, myMapF1.Extent(), Standard_False);
898       }
899
900       BoxCalculation (myMapV1, myBV1);
901       BoxCalculation (myMapE1, myBE1);
902       BoxCalculation (myMapF1, myBF1);
903
904       myIsInitS1 = Standard_True;
905     }
906
907     if (!myIsInitS2) // rebuild cached data for 2nd shape
908     {
909       if (!myMapV2.IsEmpty())
910       {
911         myBV2.Resize(1, myMapV2.Extent(), Standard_False);
912       }
913       if (!myMapE2.IsEmpty())
914       {
915         myBE2.Resize(1, myMapE2.Extent(), Standard_False);
916       }
917       if (!myMapF2.IsEmpty())
918       {
919         myBF2.Resize(1, myMapF2.Extent(), Standard_False);
920       }
921
922       BoxCalculation (myMapV2, myBV2);
923       BoxCalculation (myMapE2, myBE2);
924       BoxCalculation (myMapF2, myBF2);
925
926       myIsInitS2 = Standard_True;
927     }
928
929     if (myMapV1.Extent() && myMapV2.Extent())
930     {
931       const TopoDS_Vertex& V1 = TopoDS::Vertex(myMapV1(1));
932       const TopoDS_Vertex& V2 = TopoDS::Vertex(myMapV2(1));
933       myDistRef = DistanceInitiale(V1, V2);
934     }
935     else
936       myDistRef = 1.e30; //szv:!!!
937
938     if (!DistanceVertVert(myMapV1, myMapV2, aRootScope.Next()))
939     {
940       return Standard_False;
941     }
942     if (!DistanceMapMap(myMapV1, myMapE2, myBV1, myBE2, aRootScope.Next()))
943     {
944       return Standard_False;
945     }
946     if (!DistanceMapMap(myMapE1, myMapV2, myBE1, myBV2, aRootScope.Next()))
947     {
948       return Standard_False;
949     }
950     if (!DistanceMapMap(myMapV1, myMapF2, myBV1, myBF2, aRootScope.Next()))
951     {
952       return Standard_False;
953     }
954     if (!DistanceMapMap(myMapF1, myMapV2, myBF1, myBV2, aRootScope.Next()))
955     {
956       return Standard_False;
957     }
958     if (!DistanceMapMap(myMapE1, myMapE2, myBE1, myBE2, aRootScope.Next()))
959     {
960       return Standard_False;
961     }
962     if (!DistanceMapMap(myMapE1, myMapF2, myBE1, myBF2, aRootScope.Next()))
963     {
964       return Standard_False;
965     }
966     if (!DistanceMapMap(myMapF1, myMapE2, myBF1, myBE2, aRootScope.Next()))
967     {
968       return Standard_False;
969     }
970
971     if (Abs(myDistRef) > myEps)
972     {
973       if (!DistanceMapMap(myMapF1, myMapF2, myBF1, myBF2, aRootScope.Next()))
974       {
975         return Standard_False;
976       }
977     }
978
979     //  Modified by Sergey KHROMOV - Tue Mar  6 11:55:03 2001 Begin
980     Standard_Integer i = 1;
981     for (; i <= mySolutionsShape1.Length(); i++)
982       if (mySolutionsShape1.Value(i).Dist() > myDistRef + myEps)
983       {
984         mySolutionsShape1.Remove(i);
985         mySolutionsShape2.Remove(i);
986       }
987     //  Modified by Sergey KHROMOV - Tue Mar  6 11:55:04 2001 End
988     myIsDone = (mySolutionsShape1.Length() > 0);
989   }
990
991   return myIsDone;
992 }
993
994 //=======================================================================
995 //function : Value
996 //purpose  : 
997 //=======================================================================
998
999 Standard_Real BRepExtrema_DistShapeShape::Value() const 
1000
1001   if (!myIsDone)
1002     throw StdFail_NotDone("BRepExtrema_DistShapeShape::Value: There's no solution ");
1003
1004   return myDistRef;
1005 }
1006
1007 //=======================================================================
1008 //function : SupportOnShape1
1009 //purpose  : 
1010 //=======================================================================
1011
1012 TopoDS_Shape BRepExtrema_DistShapeShape::SupportOnShape1(const Standard_Integer N) const
1013
1014   if (!myIsDone)
1015     throw StdFail_NotDone("BRepExtrema_DistShapeShape::SupportOnShape1: There's no solution ");
1016
1017   const BRepExtrema_SolutionElem &sol = mySolutionsShape1.Value(N);
1018   switch (sol.SupportKind())
1019   {
1020     case BRepExtrema_IsVertex : return sol.Vertex();
1021     case BRepExtrema_IsOnEdge : return sol.Edge();
1022     case BRepExtrema_IsInFace : return sol.Face();
1023   }
1024   return TopoDS_Shape();
1025 }
1026
1027 //=======================================================================
1028 //function : SupportOnShape2
1029 //purpose  : 
1030 //=======================================================================
1031
1032 TopoDS_Shape BRepExtrema_DistShapeShape::SupportOnShape2(const Standard_Integer N) const 
1033
1034   if (!myIsDone)
1035     throw StdFail_NotDone("BRepExtrema_DistShapeShape::SupportOnShape2: There's no solution ");
1036
1037   const BRepExtrema_SolutionElem &sol = mySolutionsShape2.Value(N);
1038   switch (sol.SupportKind())
1039   {
1040     case BRepExtrema_IsVertex : return sol.Vertex();
1041     case BRepExtrema_IsOnEdge : return sol.Edge();
1042     case BRepExtrema_IsInFace : return sol.Face();
1043   }
1044   return TopoDS_Shape();
1045 }
1046
1047 //=======================================================================
1048 //function : ParOnEdgeS1
1049 //purpose  : 
1050 //=======================================================================
1051
1052 void BRepExtrema_DistShapeShape::ParOnEdgeS1(const Standard_Integer N, Standard_Real& t) const 
1053
1054   if (!myIsDone)
1055     throw StdFail_NotDone("BRepExtrema_DistShapeShape::ParOnEdgeS1: There's no solution");
1056
1057   const BRepExtrema_SolutionElem &sol = mySolutionsShape1.Value(N);
1058   if (sol.SupportKind() != BRepExtrema_IsOnEdge)
1059     throw BRepExtrema_UnCompatibleShape("BRepExtrema_DistShapeShape::ParOnEdgeS1: ParOnEdgeS1 is impossible without EDGE");
1060
1061   sol.EdgeParameter(t);
1062 }
1063
1064 //=======================================================================
1065 //function : ParOnEdgeS2
1066 //purpose  : 
1067 //=======================================================================
1068
1069 void BRepExtrema_DistShapeShape::ParOnEdgeS2(const Standard_Integer N,  Standard_Real& t) const 
1070
1071   if (!myIsDone)
1072     throw StdFail_NotDone("BRepExtrema_DistShapeShape::ParOnEdgeS2: There's no solution");
1073
1074   const BRepExtrema_SolutionElem &sol = mySolutionsShape2.Value(N);
1075   if (sol.SupportKind() != BRepExtrema_IsOnEdge)
1076     throw BRepExtrema_UnCompatibleShape("BRepExtrema_DistShapeShape::ParOnEdgeS2: ParOnEdgeS2 is impossible without EDGE");
1077  
1078   sol.EdgeParameter(t);
1079 }
1080
1081 //=======================================================================
1082 //function : ParOnFaceS1
1083 //purpose  : 
1084 //=======================================================================
1085
1086 void BRepExtrema_DistShapeShape::ParOnFaceS1(const Standard_Integer N,  Standard_Real& u,  Standard_Real& v) const 
1087
1088   if (!myIsDone)
1089     throw StdFail_NotDone("BRepExtrema_DistShapeShape::ParOnFaceS1: There's no solution");
1090
1091   const BRepExtrema_SolutionElem &sol = mySolutionsShape1.Value(N);
1092   if (sol.SupportKind() != BRepExtrema_IsInFace)
1093     throw BRepExtrema_UnCompatibleShape("BRepExtrema_DistShapeShape::ParOnFaceS1: ParOnFaceS1 is impossible without FACE");
1094   
1095   sol.FaceParameter(u, v);
1096 }
1097
1098 //=======================================================================
1099 //function : ParOnFaceS2
1100 //purpose  : 
1101 //=======================================================================
1102
1103 void BRepExtrema_DistShapeShape::ParOnFaceS2(const Standard_Integer N,  Standard_Real& u, Standard_Real& v) const 
1104
1105   if (!myIsDone)
1106     throw StdFail_NotDone("BRepExtrema_DistShapeShape::ParOnFaceS2: There's no solution");
1107
1108   const BRepExtrema_SolutionElem &sol = mySolutionsShape2.Value(N);
1109   if (sol.SupportKind() != BRepExtrema_IsInFace)
1110     throw BRepExtrema_UnCompatibleShape("BRepExtrema_DistShapeShape::ParOnFaceS2:ParOnFaceS2 is impossible without FACE ");
1111   
1112   sol.FaceParameter(u, v);
1113 }
1114
1115 //=======================================================================
1116 //function : Dump
1117 //purpose  : 
1118 //=======================================================================
1119
1120 void BRepExtrema_DistShapeShape::Dump(Standard_OStream& o) const 
1121 {
1122   Standard_Integer i;
1123   Standard_Real r1,r2;
1124   
1125   o<< "the distance  value is :  " << Value()<<std::endl;
1126   o<< "the number of solutions is :"<<NbSolution()<<std::endl;
1127   o<<std::endl;
1128   for (i=1;i<=NbSolution();i++)
1129   {
1130     o<<"solution number "<<i<<": "<< std::endl;
1131     o<<"the type of the solution on the first shape is " <<Standard_Integer( SupportTypeShape1(i)) <<std::endl; 
1132     o<<"the type of the solution on the second shape is "<<Standard_Integer( SupportTypeShape2(i))<< std::endl;
1133     o<< "the coordinates of  the point on the first shape are: "<<std::endl; 
1134     o<<"X=" <<PointOnShape1(i).X()<<" Y=" <<PointOnShape1(i).Y()<<" Z="<<PointOnShape1(i).Z()<<std::endl;
1135     o<< "the coordinates of  the point on the second shape are: "<<std::endl; 
1136     o<<"X="<< PointOnShape2(i).X()<< " Y="<<PointOnShape2(i).Y()<<" Z="<< PointOnShape2(i).Z()<<std::endl;
1137     
1138     switch (SupportTypeShape1(i))
1139     {
1140       case BRepExtrema_IsVertex:
1141         break;
1142       case BRepExtrema_IsOnEdge:
1143         ParOnEdgeS1(i,r1);
1144         o << "parameter on the first edge :  t= " << r1 << std::endl;
1145         break;
1146       case BRepExtrema_IsInFace:
1147         ParOnFaceS1(i,r1,r2);
1148         o << "parameters on the first face :  u= " << r1 << " v=" <<  r2 << std::endl;
1149         break;
1150     }
1151     switch (SupportTypeShape2(i))
1152     {
1153       case BRepExtrema_IsVertex:
1154         break;
1155       case BRepExtrema_IsOnEdge:
1156         ParOnEdgeS2(i,r1);
1157         o << "parameter on the second edge : t=" << r1 << std::endl;
1158         break;
1159       case BRepExtrema_IsInFace:
1160         ParOnFaceS2(i,r1,r2);
1161         o << "parameters on the second face : u= " << r1 << " v=" <<  r2 << std::endl;
1162         break;
1163     }
1164     o<<std::endl;
1165   } 
1166 }