]> SALOME platform Git repositories - modules/geom.git/blob - src/GEOMImpl/GEOMImpl_IBlocksOperations.cxx
Salome HOME
XAO test fixed to work in salome test
[modules/geom.git] / src / GEOMImpl / GEOMImpl_IBlocksOperations.cxx
1 // Copyright (C) 2007-2016  CEA/DEN, EDF R&D, OPEN CASCADE
2 //
3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
5 //
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License, or (at your option) any later version.
10 //
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14 // Lesser General Public License for more details.
15 //
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
19 //
20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
21 //
22
23 #ifdef WIN32
24 #pragma warning( disable:4786 )
25 #endif
26
27 #include <Standard_Stream.hxx>
28
29 #include <GEOMImpl_IBlocksOperations.hxx>
30
31 #include <GEOMImpl_Types.hxx>
32
33 #include <GEOMImpl_BlockDriver.hxx>
34 #include <GEOMImpl_IBlocks.hxx>
35 #include <GEOMImpl_IBlockTrsf.hxx>
36 #include <GEOMImpl_Block6Explorer.hxx>
37
38 #include <GEOMUtils.hxx>
39
40 #include <GEOM_Function.hxx>
41 #include <GEOM_PythonDump.hxx>
42
43 #include <GEOMAlgo_GlueAnalyser.hxx>
44 #include <GEOMAlgo_CoupleOfShapes.hxx>
45 #include <GEOMAlgo_ListOfCoupleOfShapes.hxx>
46 #include <GEOMAlgo_ListIteratorOfListOfCoupleOfShapes.hxx>
47 #include <BlockFix_CheckTool.hxx>
48
49 #include <Basics_OCCTVersion.hxx>
50
51 #include "utilities.h"
52 #include <OpUtil.hxx>
53 #include <Utils_ExceptHandlers.hxx>
54
55 #include <TFunction_DriverTable.hxx>
56 #include <TFunction_Driver.hxx>
57 #include <TDataStd_Integer.hxx>
58 #include <TDF_Tool.hxx>
59
60 #include <BRep_Tool.hxx>
61 #include <BRep_Builder.hxx>
62 #include <BRepTools.hxx>
63 #include <BRepTools_WireExplorer.hxx>
64 #include <BRepGProp.hxx>
65 #include <BRepBndLib.hxx>
66 #include <BRepAdaptor_Surface.hxx>
67 #include <BRepClass_FaceClassifier.hxx>
68 #include <BRepClass3d_SolidClassifier.hxx>
69 #include <BRepExtrema_DistShapeShape.hxx>
70
71 #include <TopAbs.hxx>
72 #include <TopoDS.hxx>
73 #include <TopoDS_Edge.hxx>
74 #include <TopoDS_Vertex.hxx>
75 #include <TopoDS_Compound.hxx>
76 #include <TopoDS_Iterator.hxx>
77 #include <TopExp.hxx>
78 #include <TopExp_Explorer.hxx>
79 #include <TopTools_MapOfShape.hxx>
80 #include <TopTools_Array1OfShape.hxx>
81 #include <TopTools_IndexedMapOfShape.hxx>
82 #include <TopTools_DataMapOfShapeInteger.hxx>
83 #include <TopTools_ListIteratorOfListOfShape.hxx>
84 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
85 #include <TopTools_DataMapIteratorOfDataMapOfShapeInteger.hxx>
86
87 #include <Bnd_Box.hxx>
88 #include <GProp_GProps.hxx>
89
90 #include <Geom_Curve.hxx>
91 #include <Geom_Surface.hxx>
92 #include <ShapeAnalysis_Surface.hxx>
93
94 #include <TColStd_MapOfInteger.hxx>
95 #include <TColStd_Array1OfReal.hxx>
96 #include <TColStd_Array1OfInteger.hxx>
97 #include <TColStd_Array2OfInteger.hxx>
98
99 //#include <OSD_Timer.hxx>
100
101 #include <Precision.hxx>
102
103 #include <Standard_Failure.hxx>
104 #include <Standard_ErrorHandler.hxx> // CAREFUL ! position of this file is critic : see Lucien PIGNOLONI / OCC
105
106
107 /**
108  * This function returns Standard_True if the face is quadrangular. It means
109  * that it has only 1 wire with 4 edges. If there are more then 4 edges in
110  * the wire and theToleranceC1 is not negative the new implementation is used.
111  * According to it the face is quadrangular if it is quadrangular according to
112  * an old implementation or if it has a single wire with more then 4 edges
113  * that form exactly 4 bounds of C1 continuity with the given tolerance.
114  *
115  * \param theFace the face to be checked
116  * \param theToleranceC1 if negative, it is not used; otherwise it is used
117  *        to check if two neighbor edges of face have C1 continuity.
118  * \return Standard_True if the face is quadrangular; Standard_False otherwise.
119  */
120 static Standard_Boolean IsQuadrangle(const TopoDS_Face   &theFace,
121                                      const Standard_Real  theToleranceC1)
122 {
123   TopExp_Explorer aFExp (theFace, TopAbs_WIRE);
124
125   if (!aFExp.More()) {
126     // no wire in the face
127     return Standard_False;
128   }
129
130   TopoDS_Shape aWire = aFExp.Current();
131
132   aFExp.Next();
133
134   if (aFExp.More()) {
135     // multiple wires in the face
136     return Standard_False;
137   }
138
139   // Check number of edges in the face
140   Standard_Integer    aNbEdges = 0;
141   TopTools_MapOfShape aMapEdges;
142   TopExp_Explorer     aWExp(aWire, TopAbs_EDGE);
143
144   for (; aWExp.More(); aWExp.Next()) {
145     if (aMapEdges.Add(aWExp.Current())) {
146       aNbEdges++;
147
148       if (aNbEdges > 4) {
149         break;
150       }
151     }
152   }
153
154   if (aNbEdges < 4) {
155     return Standard_False;
156   }
157
158   if (aNbEdges > 4) {
159     if (theToleranceC1 < 0.) {
160       return Standard_False;
161     }
162
163     // Check if a wire has 4 bounds of C1 continuity.
164     BRepTools_WireExplorer aWireExp(TopoDS::Wire(aWire), theFace);
165     TopTools_ListOfShape   anEdges;
166
167     for (aNbEdges = 0; aWireExp.More(); aWireExp.Next()) {
168       const TopoDS_Edge &anEdge = aWireExp.Current();
169
170       // Skip degenerated edges.
171       if (!BRep_Tool::Degenerated(anEdge)) {
172         anEdges.Append(anEdge);
173         ++aNbEdges;
174       }
175     }
176
177     if (aNbEdges < 4) {
178       return Standard_False;
179     }
180
181     // Compute number of sharp corners.
182     anEdges.Append(anEdges.First()); // To make a loop.
183
184     TopTools_ListIteratorOfListOfShape anIter(anEdges);
185     Standard_Real      aPar[2];
186     Standard_Integer   aNbCorners = 0;
187     TopoDS_Edge        anEdge1    = TopoDS::Edge(anEdges.First());
188     Handle(Geom_Curve) aCurve1    = BRep_Tool::Curve(anEdge1, aPar[0], aPar[1]);
189     Handle(Geom_Curve) aCurve2;
190     TopoDS_Edge        anEdge2;
191     TopoDS_Vertex      aCommonVtx;
192     gp_Pnt             aPnt;
193     gp_Vec             aVec1;
194     gp_Vec             aVec2;
195     Standard_Boolean   isReversed1 = (anEdge1.Orientation() == TopAbs_REVERSED);
196     Standard_Boolean   isReversed2;
197
198     for (anIter.Next(); anIter.More(); anIter.Next()) {
199       TopoDS_Edge anEdge2 = TopoDS::Edge(anIter.Value());
200
201       if (!TopExp::CommonVertex(anEdge1, anEdge2, aCommonVtx)) {
202         // NEVERREACHED
203         return Standard_False;
204       }
205
206       // Check the angle between tangent vectors of 2 curves at this point.
207       Standard_Real aParam1 = BRep_Tool::Parameter(aCommonVtx, anEdge1);
208       Standard_Real aParam2 = BRep_Tool::Parameter(aCommonVtx, anEdge2);
209
210       aCurve2     = BRep_Tool::Curve(anEdge2, aPar[0], aPar[1]);
211       isReversed2 = (anEdge2.Orientation() == TopAbs_REVERSED);
212       aCurve1->D1(aParam1, aPnt, aVec1);
213       aCurve2->D1(aParam2, aPnt, aVec2);
214
215       if (isReversed1) {
216         aVec1.Reverse();
217       }
218
219       if (isReversed2) {
220         aVec2.Reverse();
221       }
222       const Standard_Real anAngle = aVec1.Angle(aVec2);
223
224       if (anAngle > theToleranceC1) {
225         ++aNbCorners;
226
227         if (aNbCorners > 4) {
228           break;
229         }
230       }
231
232       // Go to the next couple of edges.
233       anEdge1     = anEdge2;
234       aCurve1     = aCurve2;
235       isReversed1 = isReversed2;
236     }
237
238     // Check the total number of corners.
239     if (aNbCorners != 4) {
240       return Standard_False;
241     }
242   }
243
244   return Standard_True;
245 }
246
247 //=============================================================================
248 /*!
249  *   constructor:
250  */
251 //=============================================================================
252 GEOMImpl_IBlocksOperations::GEOMImpl_IBlocksOperations (GEOM_Engine* theEngine, int theDocID)
253 : GEOM_IOperations(theEngine, theDocID)
254 {
255   MESSAGE("GEOMImpl_IBlocksOperations::GEOMImpl_IBlocksOperations");
256 }
257
258 //=============================================================================
259 /*!
260  *  destructor
261  */
262 //=============================================================================
263 GEOMImpl_IBlocksOperations::~GEOMImpl_IBlocksOperations()
264 {
265   MESSAGE("GEOMImpl_IBlocksOperations::~GEOMImpl_IBlocksOperations");
266 }
267
268
269 //=============================================================================
270 /*!
271  *  MakeQuad
272  */
273 //=============================================================================
274 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad
275                      (Handle(GEOM_Object) theEdge1, Handle(GEOM_Object) theEdge2,
276                       Handle(GEOM_Object) theEdge3, Handle(GEOM_Object) theEdge4)
277 {
278   SetErrorCode(KO);
279
280   if (theEdge1.IsNull() || theEdge2.IsNull() ||
281       theEdge3.IsNull() || theEdge4.IsNull()) return NULL;
282
283   //Add a new Face object
284   Handle(GEOM_Object) aFace = GetEngine()->AddObject(GetDocID(), GEOM_FACE);
285
286   //Add a new Face function
287   Handle(GEOM_Function) aFunction =
288     aFace->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_FACE_FOUR_EDGES);
289
290   //Check if the function is set correctly
291   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
292
293   GEOMImpl_IBlocks aPI (aFunction);
294
295   Handle(GEOM_Function) aRef1 = theEdge1->GetLastFunction();
296   Handle(GEOM_Function) aRef2 = theEdge2->GetLastFunction();
297   Handle(GEOM_Function) aRef3 = theEdge3->GetLastFunction();
298   Handle(GEOM_Function) aRef4 = theEdge4->GetLastFunction();
299   if (aRef1.IsNull() || aRef2.IsNull() ||
300       aRef3.IsNull() || aRef4.IsNull()) return NULL;
301
302   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
303   aShapesSeq->Append(aRef1);
304   aShapesSeq->Append(aRef2);
305   aShapesSeq->Append(aRef3);
306   aShapesSeq->Append(aRef4);
307
308   aPI.SetShapes(aShapesSeq);
309
310   //Compute the Face value
311   try {
312     OCC_CATCH_SIGNALS;
313     if (!GetSolver()->ComputeFunction(aFunction)) {
314       SetErrorCode("Block driver failed to compute a face");
315       return NULL;
316     }
317   }
318   catch (Standard_Failure) {
319     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
320     SetErrorCode(aFail->GetMessageString());
321     return NULL;
322   }
323
324   //Make a Python command
325   GEOM::TPythonDump(aFunction) << aFace << " = geompy.MakeQuad("
326     << theEdge1 << ", " << theEdge2 << ", " << theEdge3 << ", " << theEdge4 << ")";
327
328   SetErrorCode(OK);
329   return aFace;
330 }
331
332 //=============================================================================
333 /*!
334  *  MakeQuad2Edges
335  */
336 //=============================================================================
337 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad2Edges
338                      (Handle(GEOM_Object) theEdge1, Handle(GEOM_Object) theEdge2)
339 {
340   SetErrorCode(KO);
341
342   if (theEdge1.IsNull() || theEdge2.IsNull()) return NULL;
343
344   //Add a new Face object
345   Handle(GEOM_Object) aFace = GetEngine()->AddObject(GetDocID(), GEOM_FACE);
346
347   //Add a new Face function
348   Handle(GEOM_Function) aFunction =
349     aFace->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_FACE_TWO_EDGES);
350
351   //Check if the function is set correctly
352   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
353
354   GEOMImpl_IBlocks aPI (aFunction);
355
356   Handle(GEOM_Function) aRef1 = theEdge1->GetLastFunction();
357   Handle(GEOM_Function) aRef2 = theEdge2->GetLastFunction();
358   if (aRef1.IsNull() || aRef2.IsNull()) return NULL;
359
360   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
361   aShapesSeq->Append(aRef1);
362   aShapesSeq->Append(aRef2);
363
364   aPI.SetShapes(aShapesSeq);
365
366   //Compute the Face value
367   try {
368     OCC_CATCH_SIGNALS;
369     if (!GetSolver()->ComputeFunction(aFunction)) {
370       SetErrorCode("Block driver failed to compute a face");
371       return NULL;
372     }
373   }
374   catch (Standard_Failure) {
375     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
376     SetErrorCode(aFail->GetMessageString());
377     return NULL;
378   }
379
380   //Make a Python command
381   GEOM::TPythonDump(aFunction) << aFace << " = geompy.MakeQuad2Edges("
382                                << theEdge1 << ", " << theEdge2 << ")";
383
384   SetErrorCode(OK);
385   return aFace;
386 }
387
388 //=============================================================================
389 /*!
390  *  MakeQuad4Vertices
391  */
392 //=============================================================================
393 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeQuad4Vertices
394                      (Handle(GEOM_Object) thePnt1, Handle(GEOM_Object) thePnt2,
395                       Handle(GEOM_Object) thePnt3, Handle(GEOM_Object) thePnt4)
396 {
397   SetErrorCode(KO);
398
399   if (thePnt1.IsNull() || thePnt2.IsNull() ||
400       thePnt3.IsNull() || thePnt4.IsNull()) return NULL;
401
402   //Add a new Face object
403   Handle(GEOM_Object) aFace = GetEngine()->AddObject(GetDocID(), GEOM_FACE);
404
405   //Add a new Face function
406   Handle(GEOM_Function) aFunction =
407     aFace->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_FACE_FOUR_PNT);
408
409   //Check if the function is set correctly
410   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
411
412   GEOMImpl_IBlocks aPI (aFunction);
413
414   Handle(GEOM_Function) aRef1 = thePnt1->GetLastFunction();
415   Handle(GEOM_Function) aRef2 = thePnt2->GetLastFunction();
416   Handle(GEOM_Function) aRef3 = thePnt3->GetLastFunction();
417   Handle(GEOM_Function) aRef4 = thePnt4->GetLastFunction();
418   if (aRef1.IsNull() || aRef2.IsNull() ||
419       aRef3.IsNull() || aRef4.IsNull()) return NULL;
420
421   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
422   aShapesSeq->Append(aRef1);
423   aShapesSeq->Append(aRef2);
424   aShapesSeq->Append(aRef3);
425   aShapesSeq->Append(aRef4);
426
427   aPI.SetShapes(aShapesSeq);
428
429   //Compute the Face value
430   try {
431     OCC_CATCH_SIGNALS;
432     if (!GetSolver()->ComputeFunction(aFunction)) {
433       SetErrorCode("Block driver failed to compute a face");
434       return NULL;
435     }
436   }
437   catch (Standard_Failure) {
438     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
439     SetErrorCode(aFail->GetMessageString());
440     return NULL;
441   }
442
443   //Make a Python command
444   GEOM::TPythonDump(aFunction) << aFace << " = geompy.MakeQuad4Vertices("
445     << thePnt1 << ", " << thePnt2 << ", " << thePnt3 << ", " << thePnt4 << ")";
446
447   SetErrorCode(OK);
448   return aFace;
449 }
450
451 //=============================================================================
452 /*!
453  *  MakeHexa
454  */
455 //=============================================================================
456 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeHexa
457                      (Handle(GEOM_Object) theFace1, Handle(GEOM_Object) theFace2,
458                       Handle(GEOM_Object) theFace3, Handle(GEOM_Object) theFace4,
459                       Handle(GEOM_Object) theFace5, Handle(GEOM_Object) theFace6)
460 {
461   SetErrorCode(KO);
462
463   if (theFace1.IsNull() || theFace2.IsNull() ||
464       theFace3.IsNull() || theFace4.IsNull() ||
465       theFace5.IsNull() || theFace6.IsNull()) return NULL;
466
467   //Add a new Solid object
468   Handle(GEOM_Object) aBlock = GetEngine()->AddObject(GetDocID(), GEOM_BLOCK);
469
470   //Add a new Block function
471   Handle(GEOM_Function) aFunction =
472     aBlock->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_SIX_FACES);
473
474   //Check if the function is set correctly
475   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
476
477   GEOMImpl_IBlocks aPI (aFunction);
478
479   Handle(GEOM_Function) aRef1 = theFace1->GetLastFunction();
480   Handle(GEOM_Function) aRef2 = theFace2->GetLastFunction();
481   Handle(GEOM_Function) aRef3 = theFace3->GetLastFunction();
482   Handle(GEOM_Function) aRef4 = theFace4->GetLastFunction();
483   Handle(GEOM_Function) aRef5 = theFace5->GetLastFunction();
484   Handle(GEOM_Function) aRef6 = theFace6->GetLastFunction();
485   if (aRef1.IsNull() || aRef2.IsNull() ||
486       aRef3.IsNull() || aRef4.IsNull() ||
487       aRef5.IsNull() || aRef6.IsNull()) return NULL;
488
489   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
490   aShapesSeq->Append(aRef1);
491   aShapesSeq->Append(aRef2);
492   aShapesSeq->Append(aRef3);
493   aShapesSeq->Append(aRef4);
494   aShapesSeq->Append(aRef5);
495   aShapesSeq->Append(aRef6);
496
497   aPI.SetShapes(aShapesSeq);
498
499   //Compute the Block value
500   try {
501     OCC_CATCH_SIGNALS;
502     if (!GetSolver()->ComputeFunction(aFunction)) {
503       SetErrorCode("Block driver failed to compute a block");
504       return NULL;
505     }
506   }
507   catch (Standard_Failure) {
508     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
509     SetErrorCode(aFail->GetMessageString());
510     return NULL;
511   }
512
513   //Make a Python command
514   GEOM::TPythonDump(aFunction) << aBlock << " = geompy.MakeHexa("
515     << theFace1 << ", " << theFace2 << ", " << theFace3 << ", "
516       << theFace4 << ", " << theFace5 << ", " << theFace6 << ")";
517
518   SetErrorCode(OK);
519   return aBlock;
520 }
521
522 //=============================================================================
523 /*!
524  *  MakeHexa2Faces
525  */
526 //=============================================================================
527 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeHexa2Faces
528                    (Handle(GEOM_Object) theFace1, Handle(GEOM_Object) theFace2)
529 {
530   SetErrorCode(KO);
531
532   if (theFace1.IsNull() || theFace2.IsNull()) return NULL;
533
534   //Add a new Solid object
535   Handle(GEOM_Object) aBlock = GetEngine()->AddObject(GetDocID(), GEOM_BLOCK);
536
537   //Add a new Block function
538   Handle(GEOM_Function) aFunction =
539     aBlock->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_TWO_FACES);
540
541   //Check if the function is set correctly
542   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
543
544   GEOMImpl_IBlocks aPI (aFunction);
545
546   Handle(GEOM_Function) aRef1 = theFace1->GetLastFunction();
547   Handle(GEOM_Function) aRef2 = theFace2->GetLastFunction();
548   if (aRef1.IsNull() || aRef2.IsNull()) return NULL;
549
550   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
551   aShapesSeq->Append(aRef1);
552   aShapesSeq->Append(aRef2);
553
554   aPI.SetShapes(aShapesSeq);
555
556   //Compute the Block value
557   try {
558     OCC_CATCH_SIGNALS;
559     if (!GetSolver()->ComputeFunction(aFunction)) {
560       SetErrorCode("Block driver failed to compute a block");
561       return NULL;
562     }
563   }
564   catch (Standard_Failure) {
565     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
566     SetErrorCode(aFail->GetMessageString());
567     return NULL;
568   }
569
570   //Make a Python command
571   GEOM::TPythonDump(aFunction) << aBlock << " = geompy.MakeHexa2Faces("
572                                << theFace1 << ", " << theFace2 << ")";
573
574   SetErrorCode(OK);
575   return aBlock;
576 }
577
578 //=============================================================================
579 /*!
580  *  MakeBlockCompound
581  */
582 //=============================================================================
583 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeBlockCompound
584                                               (Handle(GEOM_Object) theCompound)
585 {
586   SetErrorCode(KO);
587
588   if (theCompound.IsNull()) return NULL;
589
590   //Add a new object
591   Handle(GEOM_Object) aBlockComp = GetEngine()->AddObject(GetDocID(), GEOM_COMPOUND);
592
593   //Add a new BlocksComp function
594   Handle(GEOM_Function) aFunction =
595     aBlockComp->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_COMPOUND_GLUE);
596
597   //Check if the function is set correctly
598   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
599
600   GEOMImpl_IBlocks aPI (aFunction);
601
602   Handle(GEOM_Function) aRef1 = theCompound->GetLastFunction();
603   if (aRef1.IsNull()) return NULL;
604
605   Handle(TColStd_HSequenceOfTransient) aShapesSeq = new TColStd_HSequenceOfTransient;
606   aShapesSeq->Append(aRef1);
607
608   aPI.SetShapes(aShapesSeq);
609
610   //Compute the Blocks Compound value
611   try {
612     OCC_CATCH_SIGNALS;
613     if (!GetSolver()->ComputeFunction(aFunction)) {
614       SetErrorCode("Block driver failed to compute a blocks compound");
615       return NULL;
616     }
617   }
618   catch (Standard_Failure) {
619     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
620     SetErrorCode(aFail->GetMessageString());
621     return NULL;
622   }
623
624   //Make a Python command
625   GEOM::TPythonDump(aFunction) << aBlockComp
626     << " = geompy.MakeBlockCompound(" << theCompound << ")";
627
628   SetErrorCode(OK);
629   return aBlockComp;
630 }
631
632 //=============================================================================
633 /*!
634  *  GetPoint
635  */
636 //=============================================================================
637 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetPoint
638                                                (Handle(GEOM_Object) theShape,
639                                                 const Standard_Real theX,
640                                                 const Standard_Real theY,
641                                                 const Standard_Real theZ,
642                                                 const Standard_Real theEpsilon)
643 {
644   SetErrorCode(KO);
645
646   //New Point object
647   Handle(GEOM_Object) aResult;
648
649   // Arguments
650   if (theShape.IsNull()) return NULL;
651
652   TopoDS_Shape aBlockOrComp = theShape->GetValue();
653   if (aBlockOrComp.IsNull()) {
654     SetErrorCode("Given shape is null");
655     return NULL;
656   }
657
658   //Compute the Vertex value
659   gp_Pnt P (theX, theY, theZ);
660   Standard_Real eps = Max(theEpsilon, Precision::Confusion());
661
662   TopoDS_Shape V;
663   Standard_Integer isFound = 0;
664   TopTools_MapOfShape mapShape;
665   TopExp_Explorer exp (aBlockOrComp, TopAbs_VERTEX);
666
667   for (; exp.More(); exp.Next()) {
668     if (mapShape.Add(exp.Current())) {
669       TopoDS_Vertex aVi = TopoDS::Vertex(exp.Current());
670       gp_Pnt aPi = BRep_Tool::Pnt(aVi);
671       if (aPi.Distance(P) < eps) {
672         V = aVi;
673         isFound++;
674       }
675     }
676   }
677
678   if (isFound == 0) {
679     SetErrorCode("Vertex has not been found");
680     return NULL;
681   } else if (isFound > 1) {
682     SetErrorCode("Multiple vertices found by the given coordinates and epsilon");
683     return NULL;
684   } else {
685     TopTools_IndexedMapOfShape anIndices;
686     TopExp::MapShapes(aBlockOrComp, anIndices);
687     Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
688     anArray->SetValue(1, anIndices.FindIndex(V));
689     aResult = GetEngine()->AddSubShape(theShape, anArray);
690   }
691
692   //The GetPoint() doesn't change object so no new function is required.
693   Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
694
695   //Make a Python command
696   GEOM::TPythonDump(aFunction, /*append=*/true)
697     << aResult << " = geompy.GetPoint(" << theShape << ", "
698     << theX << ", " << theY << ", " << theZ << ", " << theEpsilon << ")";
699
700   SetErrorCode(OK);
701   return aResult;
702 }
703
704 //=============================================================================
705 /*!
706  *  GetVertexNearPoint
707  */
708 //=============================================================================
709 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetVertexNearPoint
710                                                (Handle(GEOM_Object) theShape,
711                                                 Handle(GEOM_Object) thePoint)
712 {
713   SetErrorCode(KO);
714
715   // New Point object
716   Handle(GEOM_Object) aResult;
717
718   // Arguments
719   if (theShape.IsNull() || thePoint.IsNull()) return NULL;
720
721   TopoDS_Shape aBlockOrComp = theShape->GetValue();
722   TopoDS_Shape aPoint       = thePoint->GetValue();
723   if (aBlockOrComp.IsNull() || aPoint.IsNull()) {
724     SetErrorCode("Given shape is null");
725     return NULL;
726   }
727
728   if (aPoint.ShapeType() != TopAbs_VERTEX) {
729     SetErrorCode("Element for vertex identification is not a vertex");
730     return NULL;
731   }
732
733   TopoDS_Vertex aVert = TopoDS::Vertex(aPoint);
734   gp_Pnt aP = BRep_Tool::Pnt(aVert);
735
736   // Compute the Vertex value
737   TopoDS_Shape V;
738   bool isFound = false;
739   Standard_Real aDist = RealLast();
740   TopTools_MapOfShape mapShape;
741
742   TopExp_Explorer exp (aBlockOrComp, TopAbs_VERTEX);
743   for (; exp.More(); exp.Next()) {
744     if (mapShape.Add(exp.Current())) {
745       TopoDS_Vertex aVi = TopoDS::Vertex(exp.Current());
746       gp_Pnt aPi = BRep_Tool::Pnt(aVi);
747       Standard_Real aDisti = aPi.Distance(aP);
748       if (aDisti < aDist) {
749         V = aVi;
750         aDist = aDisti;
751         isFound = true;
752       }
753     }
754   }
755
756   if (!isFound) {
757     SetErrorCode("Vertex has not been found");
758     return NULL;
759   }
760
761   TopTools_IndexedMapOfShape anIndices;
762   TopExp::MapShapes(aBlockOrComp, anIndices);
763   Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
764   anArray->SetValue(1, anIndices.FindIndex(V));
765   aResult = GetEngine()->AddSubShape(theShape, anArray);
766
767   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
768
769   // Make a Python command
770   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetVertexNearPoint("
771                                << theShape << ", " << thePoint << ")";
772
773   SetErrorCode(OK);
774   return aResult;
775 }
776
777 //=============================================================================
778 /*!
779  *  GetEdge
780  */
781 //=============================================================================
782 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetEdge
783                                                 (Handle(GEOM_Object) theShape,
784                                                  Handle(GEOM_Object) thePoint1,
785                                                  Handle(GEOM_Object) thePoint2)
786 {
787   SetErrorCode(KO);
788
789   //New Edge object
790   Handle(GEOM_Object) aResult;
791
792   // Arguments
793   if (theShape.IsNull() || thePoint1.IsNull() || thePoint2.IsNull()) return NULL;
794
795   TopoDS_Shape aBlockOrComp = theShape->GetValue();
796   if (aBlockOrComp.IsNull()) {
797     SetErrorCode("Given shape is null");
798     return NULL;
799   }
800
801   TopoDS_Shape anArg1 = thePoint1->GetValue();
802   TopoDS_Shape anArg2 = thePoint2->GetValue();
803   if (anArg1.IsNull() || anArg2.IsNull()) {
804     SetErrorCode("Null shape is given as argument");
805     return NULL;
806   }
807   if (anArg1.ShapeType() != TopAbs_VERTEX ||
808       anArg2.ShapeType() != TopAbs_VERTEX) {
809     SetErrorCode("Element for edge identification is not a vertex");
810     return NULL;
811   }
812
813   //Compute the Edge value
814   try {
815     OCC_CATCH_SIGNALS;
816     TopTools_IndexedDataMapOfShapeListOfShape MVE;
817     GEOMImpl_Block6Explorer::MapShapesAndAncestors
818       (aBlockOrComp, TopAbs_VERTEX, TopAbs_EDGE, MVE);
819
820     TopoDS_Shape V1,V2;
821     Standard_Integer ish, ext = MVE.Extent();
822
823     if (MVE.Contains(anArg1)) {
824       V1 = anArg1;
825     } else {
826       for (ish = 1; ish <= ext; ish++) {
827         TopoDS_Shape aShi = MVE.FindKey(ish);
828         if (BRepTools::Compare(TopoDS::Vertex(anArg1), TopoDS::Vertex(aShi))) {
829           V1 = aShi;
830           break;
831         }
832       }
833     }
834
835     if (MVE.Contains(anArg2)) {
836       V2 = anArg2;
837     } else {
838       for (ish = 1; ish <= ext; ish++) {
839         TopoDS_Shape aShi = MVE.FindKey(ish);
840         if (BRepTools::Compare(TopoDS::Vertex(anArg2), TopoDS::Vertex(aShi))) {
841           V2 = aShi;
842           break;
843         }
844       }
845     }
846
847     if (V1.IsNull() || V2.IsNull()) {
848       SetErrorCode("The given vertex does not belong to the shape");
849       return NULL;
850     }
851
852     TopoDS_Shape anEdge;
853     Standard_Integer isFound =
854       GEOMImpl_Block6Explorer::FindEdge(anEdge, V1, V2, MVE, Standard_True);
855     if (isFound == 0) {
856       SetErrorCode("The given vertices do not belong to one edge of the given shape");
857       return NULL;
858     } else if (isFound > 1) {
859       SetErrorCode("Multiple edges found by the given vertices of the shape");
860       return NULL;
861     } else {
862       TopTools_IndexedMapOfShape anIndices;
863       TopExp::MapShapes(aBlockOrComp, anIndices);
864       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
865       anArray->SetValue(1, anIndices.FindIndex(anEdge));
866       aResult = GetEngine()->AddSubShape(theShape, anArray);
867     }
868   } catch (Standard_Failure) {
869     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
870     SetErrorCode(aFail->GetMessageString());
871     return NULL;
872   }
873
874   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
875
876   //Make a Python command
877   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetEdge("
878     << theShape << ", " << thePoint1 << ", " << thePoint2 << ")";
879
880   SetErrorCode(OK);
881   return aResult;
882 }
883
884 //=============================================================================
885 /*!
886  *  GetEdgeNearPoint
887  */
888 //=============================================================================
889 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetEdgeNearPoint
890                                                 (Handle(GEOM_Object) theShape,
891                                                  Handle(GEOM_Object) thePoint)
892 {
893   SetErrorCode(KO);
894
895   //New object
896   Handle(GEOM_Object) aResult;
897
898   // Arguments
899   if (theShape.IsNull() || thePoint.IsNull()) return NULL;
900
901   TopoDS_Shape aBlockOrComp = theShape->GetValue();
902   if (aBlockOrComp.IsNull()) {
903     SetErrorCode("Given shape is null");
904     return NULL;
905   }
906
907   TopoDS_Shape anArg = thePoint->GetValue();
908   if (anArg.IsNull()) {
909     SetErrorCode("Null shape is given as argument");
910     return NULL;
911   }
912   if (anArg.ShapeType() != TopAbs_VERTEX) {
913     SetErrorCode("Element for edge identification is not a vertex");
914     return NULL;
915   }
916
917   //Compute the Edge value
918   try {
919     OCC_CATCH_SIGNALS;
920     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
921     TopoDS_Shape aShape = GEOMUtils::GetEdgeNearPoint(aBlockOrComp, aVert);
922
923     TopTools_IndexedMapOfShape anIndices;
924     TopExp::MapShapes(aBlockOrComp, anIndices);
925     Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
926     anArray->SetValue(1, anIndices.FindIndex(aShape));
927     aResult = GetEngine()->AddSubShape(theShape, anArray);
928   }
929   catch (Standard_Failure) {
930     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
931     SetErrorCode(aFail->GetMessageString());
932     return NULL;
933   }
934
935   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
936
937   //Make a Python command
938   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetEdgeNearPoint("
939                                << theShape << ", " << thePoint << ")";
940
941   SetErrorCode(OK);
942   return aResult;
943 }
944
945 //=============================================================================
946 /*!
947  *  GetFaceByPoints
948  */
949 //=============================================================================
950 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByPoints
951                                                 (Handle(GEOM_Object) theShape,
952                                                  Handle(GEOM_Object) thePoint1,
953                                                  Handle(GEOM_Object) thePoint2,
954                                                  Handle(GEOM_Object) thePoint3,
955                                                  Handle(GEOM_Object) thePoint4)
956 {
957   SetErrorCode(KO);
958
959   //New object
960   Handle(GEOM_Object) aResult;
961
962   // Arguments
963   if (theShape.IsNull() ||
964       thePoint1.IsNull() || thePoint2.IsNull() ||
965       thePoint3.IsNull() || thePoint4.IsNull()) return NULL;
966
967   TopoDS_Shape aBlockOrComp = theShape->GetValue();
968   if (aBlockOrComp.IsNull()) {
969     SetErrorCode("Block or compound is null");
970     return NULL;
971   }
972
973   TopoDS_Shape anArg1 = thePoint1->GetValue();
974   TopoDS_Shape anArg2 = thePoint2->GetValue();
975   TopoDS_Shape anArg3 = thePoint3->GetValue();
976   TopoDS_Shape anArg4 = thePoint4->GetValue();
977   if (anArg1.IsNull() || anArg2.IsNull() ||
978       anArg3.IsNull() || anArg4.IsNull()) {
979     SetErrorCode("Null shape is given as argument");
980     return NULL;
981   }
982   if (anArg1.ShapeType() != TopAbs_VERTEX ||
983       anArg2.ShapeType() != TopAbs_VERTEX ||
984       anArg3.ShapeType() != TopAbs_VERTEX ||
985       anArg4.ShapeType() != TopAbs_VERTEX) {
986     SetErrorCode("Element for face identification is not a vertex");
987     return NULL;
988   }
989
990   //Compute the Face value
991   try {
992     OCC_CATCH_SIGNALS;
993     TopoDS_Shape aShape;
994
995     TopTools_IndexedDataMapOfShapeListOfShape MVF;
996     GEOMImpl_Block6Explorer::MapShapesAndAncestors(aBlockOrComp, TopAbs_VERTEX, TopAbs_FACE, MVF);
997
998     TopoDS_Shape V1,V2,V3,V4;
999     Standard_Integer ish, ext = MVF.Extent();
1000
1001     if (MVF.Contains(anArg1)) {
1002       V1 = anArg1;
1003     } else {
1004       for (ish = 1; ish <= ext; ish++) {
1005         TopoDS_Shape aShi = MVF.FindKey(ish);
1006         if (BRepTools::Compare(TopoDS::Vertex(anArg1), TopoDS::Vertex(aShi))) {
1007           V1 = aShi;
1008           break;
1009         }
1010       }
1011     }
1012
1013     if (MVF.Contains(anArg2)) {
1014       V2 = anArg2;
1015     } else {
1016       for (ish = 1; ish <= ext; ish++) {
1017         TopoDS_Shape aShi = MVF.FindKey(ish);
1018         if (BRepTools::Compare(TopoDS::Vertex(anArg2), TopoDS::Vertex(aShi))) {
1019           V2 = aShi;
1020           break;
1021         }
1022       }
1023     }
1024
1025     if (MVF.Contains(anArg3)) {
1026       V3 = anArg3;
1027     } else {
1028       for (ish = 1; ish <= ext; ish++) {
1029         TopoDS_Shape aShi = MVF.FindKey(ish);
1030         if (BRepTools::Compare(TopoDS::Vertex(anArg3), TopoDS::Vertex(aShi))) {
1031           V3 = aShi;
1032           break;
1033         }
1034       }
1035     }
1036
1037     if (MVF.Contains(anArg4)) {
1038       V4 = anArg4;
1039     } else {
1040       for (ish = 1; ish <= ext; ish++) {
1041         TopoDS_Shape aShi = MVF.FindKey(ish);
1042         if (BRepTools::Compare(TopoDS::Vertex(anArg4), TopoDS::Vertex(aShi))) {
1043           V4 = aShi;
1044           break;
1045         }
1046       }
1047     }
1048
1049     if (V1.IsNull() || V2.IsNull() || V3.IsNull() || V4.IsNull()) {
1050       SetErrorCode("The given vertex does not belong to the shape");
1051       return NULL;
1052     }
1053
1054     Standard_Integer isFound =
1055       GEOMImpl_Block6Explorer::FindFace(aShape, V1, V2, V3, V4, MVF, Standard_True);
1056     if (isFound == 0) {
1057       SetErrorCode("The given vertices do not belong to one face of the given shape");
1058       return NULL;
1059     } else if (isFound > 1) {
1060       SetErrorCode("The given vertices belong to several faces of the given shape");
1061       return NULL;
1062     } else {
1063       TopTools_IndexedMapOfShape anIndices;
1064       TopExp::MapShapes(aBlockOrComp, anIndices);
1065       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
1066       anArray->SetValue(1, anIndices.FindIndex(aShape));
1067       aResult = GetEngine()->AddSubShape(theShape, anArray);
1068     }
1069   }
1070   catch (Standard_Failure) {
1071     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1072     SetErrorCode(aFail->GetMessageString());
1073     return NULL;
1074   }
1075
1076   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1077
1078   //Make a Python command
1079   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceByPoints("
1080     << theShape << ", " << thePoint1 << ", " << thePoint2
1081       << ", " << thePoint3 << ", " << thePoint4 << ")";
1082
1083   SetErrorCode(OK);
1084   return aResult;
1085 }
1086
1087 //=============================================================================
1088 /*!
1089  *  GetFaceByEdges
1090  */
1091 //=============================================================================
1092 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByEdges
1093                                                 (Handle(GEOM_Object) theShape,
1094                                                  Handle(GEOM_Object) theEdge1,
1095                                                  Handle(GEOM_Object) theEdge2)
1096 {
1097   SetErrorCode(KO);
1098
1099   //New object
1100   Handle(GEOM_Object) aResult;
1101
1102   // Arguments
1103   if (theShape.IsNull() || theEdge1.IsNull() || theEdge2.IsNull()) return NULL;
1104
1105   TopoDS_Shape aBlockOrComp = theShape->GetValue();
1106   if (aBlockOrComp.IsNull()) {
1107     SetErrorCode("Block or compound is null");
1108     return NULL;
1109   }
1110
1111   TopoDS_Shape anArg1 = theEdge1->GetValue();
1112   TopoDS_Shape anArg2 = theEdge2->GetValue();
1113   if (anArg1.IsNull() || anArg2.IsNull()) {
1114     SetErrorCode("Null shape is given as argument");
1115     return NULL;
1116   }
1117   if (anArg1.ShapeType() != TopAbs_EDGE ||
1118       anArg2.ShapeType() != TopAbs_EDGE) {
1119     SetErrorCode("Element for face identification is not an edge");
1120     return NULL;
1121   }
1122
1123   //Compute the Face value
1124   try {
1125     OCC_CATCH_SIGNALS;
1126     TopoDS_Shape aShape;
1127
1128     TopTools_IndexedDataMapOfShapeListOfShape MEF;
1129     GEOMImpl_Block6Explorer::MapShapesAndAncestors(aBlockOrComp, TopAbs_EDGE, TopAbs_FACE, MEF);
1130
1131     TopoDS_Shape E1,E2;
1132     Standard_Integer ish, ext = MEF.Extent();
1133
1134     if (MEF.Contains(anArg1)) {
1135       E1 = anArg1;
1136     } else {
1137       for (ish = 1; ish <= ext; ish++) {
1138         TopoDS_Shape aShi = MEF.FindKey(ish);
1139         if (GEOMImpl_Block6Explorer::IsSimilarEdges(anArg1, aShi)) {
1140           E1 = aShi;
1141         }
1142       }
1143     }
1144
1145     if (MEF.Contains(anArg2)) {
1146       E2 = anArg2;
1147     } else {
1148       for (ish = 1; ish <= ext; ish++) {
1149         TopoDS_Shape aShi = MEF.FindKey(ish);
1150         if (GEOMImpl_Block6Explorer::IsSimilarEdges(anArg2, aShi)) {
1151           E2 = aShi;
1152         }
1153       }
1154     }
1155
1156     if (E1.IsNull() || E2.IsNull()) {
1157       SetErrorCode("The given edge does not belong to the shape");
1158       return NULL;
1159     }
1160
1161     const TopTools_ListOfShape& aFacesOfE1 = MEF.FindFromKey(E1);
1162     const TopTools_ListOfShape& aFacesOfE2 = MEF.FindFromKey(E2);
1163
1164     Standard_Integer isFound = 0;
1165     TopTools_ListIteratorOfListOfShape anIterF1 (aFacesOfE1);
1166     for (; anIterF1.More(); anIterF1.Next()) {
1167
1168       TopTools_ListIteratorOfListOfShape anIterF2 (aFacesOfE2);
1169       for (; anIterF2.More(); anIterF2.Next()) {
1170
1171         if (anIterF1.Value().IsSame(anIterF2.Value())) {
1172           isFound++;
1173
1174           // Store the face, defined by two edges
1175           aShape = anIterF1.Value();
1176         }
1177       }
1178     }
1179     if (isFound == 0) {
1180       SetErrorCode("The given edges do not belong to one face of the given shape");
1181       return NULL;
1182     } else if (isFound > 1) {
1183       SetErrorCode("The given edges belong to several faces of the given shape");
1184       return NULL;
1185     } else {
1186       TopTools_IndexedMapOfShape anIndices;
1187       TopExp::MapShapes(aBlockOrComp, anIndices);
1188       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
1189       anArray->SetValue(1, anIndices.FindIndex(aShape));
1190       aResult = GetEngine()->AddSubShape(theShape, anArray);
1191     }
1192   }
1193   catch (Standard_Failure) {
1194     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1195     SetErrorCode(aFail->GetMessageString());
1196     return NULL;
1197   }
1198
1199   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1200
1201   //Make a Python command
1202   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceByEdges("
1203     << theShape << ", " << theEdge1 << ", " << theEdge2 << ")";
1204
1205   SetErrorCode(OK);
1206   return aResult;
1207 }
1208
1209 //=============================================================================
1210 /*!
1211  *  GetOppositeFace
1212  */
1213 //=============================================================================
1214 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetOppositeFace
1215                                                 (Handle(GEOM_Object) theShape,
1216                                                  Handle(GEOM_Object) theFace)
1217 {
1218   SetErrorCode(KO);
1219
1220   //New object
1221   Handle(GEOM_Object) aResult;
1222
1223   // Arguments
1224   if (theShape.IsNull() || theFace.IsNull()) return NULL;
1225
1226   TopoDS_Shape aBlockOrComp = theShape->GetValue();
1227   if (aBlockOrComp.IsNull()) {
1228     SetErrorCode("Block is null");
1229     return NULL;
1230   }
1231   if (aBlockOrComp.ShapeType() != TopAbs_SOLID) {
1232     SetErrorCode("Shape is not a block");
1233     return NULL;
1234   }
1235
1236   TopoDS_Shape anArg = theFace->GetValue();
1237   if (anArg.IsNull()) {
1238     SetErrorCode("Null shape is given as argument");
1239     return NULL;
1240   }
1241   if (anArg.ShapeType() != TopAbs_FACE) {
1242     SetErrorCode("Element for face identification is not a face");
1243     return NULL;
1244   }
1245
1246   //Compute the Face value
1247   try {
1248     OCC_CATCH_SIGNALS;
1249     TopoDS_Shape aShape;
1250
1251     GEOMImpl_Block6Explorer aBlockTool;
1252     aBlockTool.InitByBlockAndFace(aBlockOrComp, anArg);
1253     aShape = aBlockTool.GetFace(2);
1254
1255     TopTools_IndexedMapOfShape anIndices;
1256     TopExp::MapShapes(aBlockOrComp, anIndices);
1257     Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
1258     anArray->SetValue(1, anIndices.FindIndex(aShape));
1259     aResult = GetEngine()->AddSubShape(theShape, anArray);
1260   }
1261   catch (Standard_Failure) {
1262     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1263     SetErrorCode(aFail->GetMessageString());
1264     return NULL;
1265   }
1266
1267   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1268
1269   //Make a Python command
1270   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetOppositeFace("
1271                                << theShape << ", " << theFace << ")";
1272
1273   SetErrorCode(OK);
1274   return aResult;
1275 }
1276
1277 //=============================================================================
1278 /*!
1279  *  GetFaceNearPoint
1280  */
1281 //=============================================================================
1282 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceNearPoint
1283                                                 (Handle(GEOM_Object) theShape,
1284                                                  Handle(GEOM_Object) thePoint)
1285 {
1286   SetErrorCode(KO);
1287
1288   //New object
1289   Handle(GEOM_Object) aResult;
1290
1291   // Arguments
1292   if (theShape.IsNull() || thePoint.IsNull()) return NULL;
1293
1294   TopoDS_Shape aBlockOrComp = theShape->GetValue();
1295   if (aBlockOrComp.IsNull()) {
1296     SetErrorCode("Block or compound is null");
1297     return NULL;
1298   }
1299
1300   TopoDS_Shape anArg = thePoint->GetValue();
1301   if (anArg.IsNull()) {
1302     SetErrorCode("Null shape is given as argument");
1303     return NULL;
1304   }
1305   if (anArg.ShapeType() != TopAbs_VERTEX) {
1306     SetErrorCode("Element for face identification is not a vertex");
1307     return NULL;
1308   }
1309
1310   //Compute the Face value
1311   try {
1312     OCC_CATCH_SIGNALS;
1313     TopoDS_Shape aShape;
1314
1315     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
1316     gp_Pnt aPnt = BRep_Tool::Pnt(aVert);
1317     Standard_Real PX, PY, PZ;
1318     aPnt.Coord(PX, PY, PZ);
1319
1320     // 1. Classify the point relatively each face
1321     Standard_Integer nearest = 2, nbFound = 0;
1322     TopTools_DataMapOfShapeInteger mapShapeDist;
1323     TopExp_Explorer exp (aBlockOrComp, TopAbs_FACE);
1324     for (; exp.More(); exp.Next()) {
1325       TopoDS_Shape aFace = exp.Current();
1326
1327       if (!mapShapeDist.IsBound(aFace)) {
1328         Standard_Integer aDistance = 2;
1329
1330         // 1.a. Classify relatively Surface
1331         Handle(Geom_Surface) aSurf = BRep_Tool::Surface(TopoDS::Face(aFace));
1332         Handle(ShapeAnalysis_Surface) aSurfAna = new ShapeAnalysis_Surface (aSurf);
1333         gp_Pnt2d p2dOnSurf = aSurfAna->ValueOfUV(aPnt, Precision::Confusion());
1334         gp_Pnt p3dOnSurf = aSurfAna->Value(p2dOnSurf);
1335         Standard_Real aDist = p3dOnSurf.Distance(aPnt);
1336         if (aDist > Precision::Confusion()) {
1337           // OUT of Surface
1338           aDistance = 1;
1339         } else {
1340           // 1.b. Classify relatively the face itself
1341           BRepClass_FaceClassifier FC (TopoDS::Face(aFace), p2dOnSurf, Precision::Confusion());
1342           if (FC.State() == TopAbs_IN) {
1343             aDistance = -1;
1344           } else if (FC.State() == TopAbs_ON) {
1345             aDistance = 0;
1346           } else { // OUT
1347             aDistance = 1;
1348           }
1349         }
1350
1351         if (aDistance < nearest) {
1352           nearest = aDistance;
1353           aShape = aFace;
1354           nbFound = 1;
1355
1356           // A first found face, containing the point inside, will be returned.
1357           // It is the solution, if there are no
1358           // coincident or intersecting faces in the compound.
1359           if (nearest == -1) break;
1360
1361         } else if (aDistance == nearest) {
1362           nbFound++;
1363         } else {
1364         }
1365
1366         mapShapeDist.Bind(aFace, aDistance);
1367       } // if (!mapShapeDist.IsBound(aFace))
1368     }
1369
1370     // 2. Define face, containing the point or having minimum distance to it
1371     if (nbFound > 1) {
1372       if (nearest == 0) {
1373         // The point is on boundary of some faces and there are
1374         // no faces, having the point inside
1375         SetErrorCode("Multiple faces near the given point are found");
1376         return NULL;
1377
1378       } else if (nearest == 1) {
1379         // The point is outside some faces and there are
1380         // no faces, having the point inside or on boundary.
1381         // We will get a nearest face
1382         Standard_Real bigReal = RealLast();
1383         Standard_Real minDist = bigReal;
1384         TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
1385         for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
1386           if (mapShapeDistIter.Value() == 1) {
1387             TopoDS_Shape aFace = mapShapeDistIter.Key();
1388             Standard_Real aDist = bigReal;
1389
1390             // 2.a. Fast check of distance - if point projection on surface is on face
1391             Handle(Geom_Surface) aSurf = BRep_Tool::Surface(TopoDS::Face(aFace));
1392             Handle(ShapeAnalysis_Surface) aSurfAna = new ShapeAnalysis_Surface (aSurf);
1393             gp_Pnt2d p2dOnSurf = aSurfAna->ValueOfUV(aPnt, Precision::Confusion());
1394             gp_Pnt p3dOnSurf = aSurfAna->Value(p2dOnSurf);
1395             aDist = p3dOnSurf.Distance(aPnt);
1396
1397             BRepClass_FaceClassifier FC (TopoDS::Face(aFace), p2dOnSurf, Precision::Confusion());
1398             if (FC.State() == TopAbs_OUT) {
1399               if (aDist < minDist) {
1400                 // 2.b. Slow check - if point projection on surface is outside of face
1401                 BRepExtrema_DistShapeShape aDistTool (aVert, aFace);
1402                 if (!aDistTool.IsDone()) {
1403                   SetErrorCode("Can not find a distance from the given point to one of faces");
1404                   return NULL;
1405                 }
1406                 aDist = aDistTool.Value();
1407               } else {
1408                 aDist = bigReal;
1409               }
1410             }
1411
1412             if (aDist < minDist) {
1413               minDist = aDist;
1414               aShape = aFace;
1415             }
1416           }
1417         }
1418       } else { // nearest == -1
1419 //        // The point is inside some faces.
1420 //        // We will get a face with nearest center
1421 //        Standard_Real minDist = RealLast();
1422 //        TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
1423 //        for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
1424 //          if (mapShapeDistIter.Value() == -1) {
1425 //            TopoDS_Shape aFace = mapShapeDistIter.Key();
1426 //            GProp_GProps aSystem;
1427 //            BRepGProp::SurfaceProperties(aFace, aSystem);
1428 //            gp_Pnt aCenterMass = aSystem.CentreOfMass();
1429 //
1430 //            Standard_Real aDist = aCenterMass.Distance(aPnt);
1431 //            if (aDist < minDist) {
1432 //              minDist = aDist;
1433 //              aShape = aFace;
1434 //            }
1435 //          }
1436 //        }
1437       }
1438     } // if (nbFound > 1)
1439
1440     if (nbFound == 0) {
1441       SetErrorCode("There are no faces near the given point");
1442       return NULL;
1443     } else {
1444       TopTools_IndexedMapOfShape anIndices;
1445       TopExp::MapShapes(aBlockOrComp, anIndices);
1446       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
1447       anArray->SetValue(1, anIndices.FindIndex(aShape));
1448       aResult = GetEngine()->AddSubShape(theShape, anArray);
1449     }
1450   }
1451   catch (Standard_Failure) {
1452     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1453     SetErrorCode(aFail->GetMessageString());
1454     return NULL;
1455   }
1456
1457   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1458
1459   //Make a Python command
1460   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceNearPoint("
1461                                << theShape << ", " << thePoint << ")";
1462
1463   SetErrorCode(OK);
1464   return aResult;
1465 }
1466
1467 //=============================================================================
1468 /*!
1469  *  GetFaceByNormale
1470  */
1471 //=============================================================================
1472 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetFaceByNormale
1473                                                 (Handle(GEOM_Object) theShape,
1474                                                  Handle(GEOM_Object) theVector)
1475 {
1476   SetErrorCode(KO);
1477
1478   //New object
1479   Handle(GEOM_Object) aResult;
1480
1481   // Arguments
1482   if (theShape.IsNull() || theVector.IsNull()) return NULL;
1483
1484   TopoDS_Shape aBlockOrComp = theShape->GetValue();
1485   if (aBlockOrComp.IsNull()) {
1486     SetErrorCode("Block or compound is null");
1487     return NULL;
1488   }
1489
1490   TopoDS_Shape anArg = theVector->GetValue();
1491   if (anArg.IsNull()) {
1492     SetErrorCode("Null shape is given as argument");
1493     return NULL;
1494   }
1495   if (anArg.ShapeType() != TopAbs_EDGE) {
1496     SetErrorCode("Element for normale identification is not an edge");
1497     return NULL;
1498   }
1499
1500   //Compute the Face value
1501   try {
1502     OCC_CATCH_SIGNALS;
1503     TopoDS_Shape aShape;
1504
1505     TopoDS_Edge anEdge = TopoDS::Edge(anArg);
1506     TopoDS_Vertex V1, V2;
1507     TopExp::Vertices(anEdge, V1, V2, Standard_True);
1508     gp_Pnt P1 = BRep_Tool::Pnt(V1);
1509     gp_Pnt P2 = BRep_Tool::Pnt(V2);
1510     gp_Vec aVec (P1, P2);
1511     if (aVec.Magnitude() < Precision::Confusion()) {
1512       SetErrorCode("Vector with null magnitude is given");
1513       return NULL;
1514     }
1515
1516     Standard_Real minAngle = RealLast();
1517     TopTools_MapOfShape mapShape;
1518     TopExp_Explorer exp (aBlockOrComp, TopAbs_FACE);
1519     for (; exp.More(); exp.Next()) {
1520       if (mapShape.Add(exp.Current())) {
1521         TopoDS_Face aFace = TopoDS::Face(exp.Current());
1522         BRepAdaptor_Surface SF (aFace);
1523
1524         Standard_Real u, v, x;
1525
1526         // find a point on the surface to get normal direction in
1527         u = SF.FirstUParameter();
1528         x = SF.LastUParameter();
1529         if (Precision::IsInfinite(u)) {
1530           u =  (Precision::IsInfinite(x)) ? 0. : x;
1531         } else if (!Precision::IsInfinite(x)) {
1532           u = (u+x) / 2.;
1533         }
1534
1535         v = SF.FirstVParameter();
1536         x = SF.LastVParameter();
1537         if (Precision::IsInfinite(v)) {
1538           v =  (Precision::IsInfinite(x)) ? 0. : x;
1539         } else if (!Precision::IsInfinite(x)) {
1540           v = (v+x) / 2.;
1541         }
1542
1543         // compute the normal direction
1544         gp_Vec Vec1,Vec2;
1545         SF.D1(u,v,P1,Vec1,Vec2);
1546         gp_Vec V = Vec1.Crossed(Vec2);
1547         x = V.Magnitude();
1548         if (V.Magnitude() < Precision::Confusion()) {
1549           SetErrorCode("Normal vector of a face has null magnitude");
1550           return NULL;
1551         }
1552
1553         // consider the face orientation
1554         if (aFace.Orientation() == TopAbs_REVERSED ||
1555             aFace.Orientation() == TopAbs_INTERNAL) {
1556           V = - V;
1557         }
1558
1559         // compute the angle and compare with the minimal one
1560         Standard_Real anAngle = aVec.Angle(V);
1561         if (anAngle < minAngle) {
1562           minAngle = anAngle;
1563           aShape = aFace;
1564         }
1565       }
1566     }
1567
1568     if (aShape.IsNull()) {
1569       SetErrorCode("Failed to find a face by the given normale");
1570       return NULL;
1571     } else {
1572       TopTools_IndexedMapOfShape anIndices;
1573       TopExp::MapShapes(aBlockOrComp, anIndices);
1574       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
1575       anArray->SetValue(1, anIndices.FindIndex(aShape));
1576       aResult = GetEngine()->AddSubShape(theShape, anArray);
1577     }
1578   }
1579   catch (Standard_Failure) {
1580     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1581     SetErrorCode(aFail->GetMessageString());
1582     return NULL;
1583   }
1584
1585   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1586
1587   //Make a Python command
1588   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetFaceByNormale("
1589     << theShape << ", " << theVector << ")";
1590
1591   SetErrorCode(OK);
1592   return aResult;
1593 }
1594
1595 //=============================================================================
1596 /*!
1597  *  GetShapesNearPoint
1598  */
1599 //=============================================================================
1600 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetShapesNearPoint
1601                                          (Handle(GEOM_Object)    theShape,
1602                                           Handle(GEOM_Object)    thePoint,
1603                                           const Standard_Integer theShapeType,
1604                                           const Standard_Real    theConstTolerance)
1605 {
1606   SetErrorCode(KO);
1607
1608   // New object
1609   Handle(GEOM_Object) aResult;
1610
1611   // Arguments
1612   if (theShape.IsNull() || thePoint.IsNull()) return NULL;
1613
1614   TopoDS_Shape aBlockOrComp = theShape->GetValue();
1615   if (aBlockOrComp.IsNull()) {
1616     SetErrorCode("Block or compound is null");
1617     return NULL;
1618   }
1619
1620   TopoDS_Shape anArg = thePoint->GetValue();
1621   if (anArg.IsNull()) {
1622     SetErrorCode("Null shape is given as argument");
1623     return NULL;
1624   }
1625   if (anArg.ShapeType() != TopAbs_VERTEX) {
1626     SetErrorCode("Element for face identification is not a vertex");
1627     return NULL;
1628   }
1629
1630   if (theShapeType < TopAbs_SOLID || TopAbs_VERTEX < theShapeType) {
1631     SetErrorCode("Invalid type of result is requested");
1632     return NULL;
1633   }
1634   
1635   Standard_Real theTolerance = theConstTolerance;
1636   if (theTolerance < Precision::Confusion()) {
1637     theTolerance = Precision::Confusion();
1638   }
1639
1640   // Compute the result
1641   try {
1642     OCC_CATCH_SIGNALS;
1643     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
1644
1645     TopTools_MapOfShape mapShape;
1646     Standard_Integer nbEdges = 0;
1647     TopExp_Explorer exp (aBlockOrComp, TopAbs_ShapeEnum(theShapeType));
1648     for (; exp.More(); exp.Next()) {
1649       if (mapShape.Add(exp.Current())) {
1650         nbEdges++;
1651       }
1652     }
1653
1654     if (nbEdges == 0) {
1655       SetErrorCode("Given shape contains no sub-shapes of requested type");
1656       return NULL;
1657     }
1658
1659     // Calculate distances and find min
1660     mapShape.Clear();
1661     Standard_Integer ind = 1;
1662     Standard_Real aMinDist = RealLast();
1663     TopTools_Array1OfShape anEdges (1, nbEdges);
1664     TColStd_Array1OfReal aDistances (1, nbEdges);
1665     for (exp.Init(aBlockOrComp, TopAbs_ShapeEnum(theShapeType)); exp.More(); exp.Next()) {
1666       if (mapShape.Add(exp.Current())) {
1667         TopoDS_Shape anEdge = exp.Current();
1668         anEdges(ind) = anEdge;
1669
1670         BRepExtrema_DistShapeShape aDistTool (aVert, anEdges(ind));
1671         if (!aDistTool.IsDone()) {
1672           SetErrorCode("Can not find a distance from the given point to one of sub-shapes");
1673           return NULL;
1674         }
1675         aDistances(ind) = aDistTool.Value();
1676         if (aDistances(ind) < aMinDist) {
1677           aMinDist = aDistances(ind);
1678         }
1679         ind++;
1680       }
1681     }
1682
1683     if (aMinDist < RealLast()) {
1684       // Collect sub-shapes with distance < (aMinDist + theTolerance)
1685       int nbSubShapes = 0;
1686       TopTools_Array1OfShape aNearShapes (1, nbEdges);
1687       for (ind = 1; ind <= nbEdges; ind++) {
1688         if (aDistances(ind) < aMinDist + theTolerance) {
1689           nbSubShapes++;
1690           aNearShapes(nbSubShapes) = anEdges(ind);
1691         }
1692       }
1693
1694       // Add sub-shape
1695       TopTools_IndexedMapOfShape anIndices;
1696       TopExp::MapShapes(aBlockOrComp, anIndices);
1697       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger (1, nbSubShapes);
1698       for (ind = 1; ind <= nbSubShapes; ind++) {
1699         anArray->SetValue(ind, anIndices.FindIndex(aNearShapes(ind)));
1700       }
1701       aResult = GetEngine()->AddSubShape(theShape, anArray);
1702     }
1703   }
1704   catch (Standard_Failure) {
1705     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1706     SetErrorCode(aFail->GetMessageString());
1707     return NULL;
1708   }
1709
1710   if (aResult.IsNull())
1711     return NULL;
1712
1713   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
1714
1715   //Make a Python command
1716   GEOM::TPythonDump(aFunction)
1717     << aResult << " = geompy.GetShapesNearPoint(" << theShape << ", " << thePoint
1718     << ", " << TopAbs_ShapeEnum(theShapeType) << ", " << theTolerance << ")";
1719
1720   SetErrorCode(OK);
1721   return aResult;
1722 }
1723
1724 //=============================================================================
1725 /*!
1726  *  IsCompoundOfBlocks
1727  */
1728 //=============================================================================
1729 Standard_Boolean GEOMImpl_IBlocksOperations::IsCompoundOfBlocks
1730                                                 (Handle(GEOM_Object)    theCompound,
1731                                                  const Standard_Integer theMinNbFaces,
1732                                                  const Standard_Integer theMaxNbFaces,
1733                                                  Standard_Integer&      theNbBlocks)
1734 {
1735   SetErrorCode(KO);
1736   Standard_Boolean isCompOfBlocks = Standard_False;
1737   theNbBlocks = 0;
1738
1739   if (theCompound.IsNull()) return isCompOfBlocks;
1740   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
1741
1742   //Check
1743   isCompOfBlocks = Standard_True;
1744   try {
1745     OCC_CATCH_SIGNALS;
1746     TopTools_MapOfShape mapShape;
1747     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
1748     for (; exp.More(); exp.Next()) {
1749       if (mapShape.Add(exp.Current())) {
1750         TopoDS_Shape aSolid = exp.Current();
1751
1752         TopTools_MapOfShape mapFaces;
1753         TopExp_Explorer expF (aSolid, TopAbs_FACE);
1754         Standard_Integer nbFaces = 0;
1755         for (; expF.More(); expF.Next()) {
1756           if (mapFaces.Add(expF.Current())) {
1757             nbFaces++;
1758             if (nbFaces > theMaxNbFaces) {
1759               isCompOfBlocks = Standard_False;
1760               break;
1761             }
1762           }
1763         }
1764         if (nbFaces < theMinNbFaces || theMaxNbFaces < nbFaces) {
1765           isCompOfBlocks = Standard_False;
1766         } else {
1767           theNbBlocks++;
1768         }
1769       }
1770     }
1771   }
1772   catch (Standard_Failure) {
1773     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1774     SetErrorCode(aFail->GetMessageString());
1775     return isCompOfBlocks;
1776   }
1777
1778   SetErrorCode(OK);
1779   return isCompOfBlocks;
1780 }
1781
1782 //=============================================================================
1783 /*!
1784  *  Set of functions, used by CheckCompoundOfBlocks() method
1785  */
1786 //=============================================================================
1787 void GEOMImpl_IBlocksOperations::AddBlocksFrom (const TopoDS_Shape&   theShape,
1788                                                 TopTools_ListOfShape& BLO,
1789                                                 TopTools_ListOfShape& NOT,
1790                                                 TopTools_ListOfShape& EXT,
1791                                                 TopTools_ListOfShape& NOQ,
1792                                                 const Standard_Real   theToleranceC1)
1793 {
1794   TopAbs_ShapeEnum aType = theShape.ShapeType();
1795   switch (aType) {
1796   case TopAbs_COMPOUND:
1797   case TopAbs_COMPSOLID:
1798     {
1799       TopoDS_Iterator It (theShape);
1800       for (; It.More(); It.Next()) {
1801         AddBlocksFrom(It.Value(), BLO, NOT, EXT, NOQ, theToleranceC1);
1802       }
1803     }
1804     break;
1805   case TopAbs_SOLID:
1806     {
1807       // Check, if there are seam or degenerated edges
1808       BlockFix_CheckTool aTool;
1809       aTool.SetShape(theShape);
1810       aTool.SetAngTolerance(theToleranceC1);
1811       aTool.Perform();
1812       if (aTool.NbPossibleBlocks() > 0) {
1813         EXT.Append(theShape);
1814       } else {
1815         // Count faces and edges in each face to recognize blocks
1816         TopTools_MapOfShape mapFaces;
1817         Standard_Integer nbFaces = 0;
1818         Standard_Boolean hasNonQuadr = Standard_False;
1819         TopExp_Explorer expF (theShape, TopAbs_FACE);
1820
1821         for (; expF.More(); expF.Next()) {
1822           TopoDS_Face aF = TopoDS::Face(expF.Current());
1823
1824           if (mapFaces.Add(aF)) {
1825             nbFaces++;
1826
1827             if (!IsQuadrangle(aF, theToleranceC1)) {
1828               hasNonQuadr = Standard_True;
1829               NOQ.Append(aF);//0021483
1830             }
1831           }
1832         }
1833
1834         if (nbFaces == 6 && !hasNonQuadr) {
1835           BLO.Append(theShape);
1836         } else {
1837           NOT.Append(theShape);
1838         }
1839       }
1840     }
1841     break;
1842   case TopAbs_SHELL: //0021483
1843   case TopAbs_FACE: //0021483
1844     {
1845       // Count edges in each face
1846       TopTools_MapOfShape mapFaces;
1847       TopExp_Explorer expF (theShape, TopAbs_FACE);
1848       for (; expF.More(); expF.Next()) {
1849         TopoDS_Face aF = TopoDS::Face(expF.Current());
1850
1851         if (mapFaces.Add(aF)) {
1852           if (!IsQuadrangle(aF, theToleranceC1)) {
1853             NOQ.Append(aF);//0021483
1854           }
1855         }
1856       }
1857     }
1858     break;
1859   default:
1860     NOT.Append(theShape);
1861   }
1862 }
1863
1864 #define REL_NOT_CONNECTED 0
1865 #define REL_OK            1
1866 #define REL_NOT_GLUED     2
1867 #define REL_COLLISION_VV  3
1868 #define REL_COLLISION_FF  4
1869 #define REL_COLLISION_EE  5
1870 #define REL_UNKNOWN       6
1871
1872 Standard_Integer BlocksRelation (const TopoDS_Shape& theBlock1,
1873                                  const TopoDS_Shape& theBlock2)
1874 {
1875   // Compare bounding boxes before calling BRepExtrema_DistShapeShape
1876   Standard_Real Xmin1, Ymin1, Zmin1, Xmax1, Ymax1, Zmax1;
1877   Standard_Real Xmin2, Ymin2, Zmin2, Xmax2, Ymax2, Zmax2;
1878   Bnd_Box B1, B2;
1879   BRepBndLib::Add(theBlock1, B1);
1880   BRepBndLib::Add(theBlock2, B2);
1881   B1.Get(Xmin1, Ymin1, Zmin1, Xmax1, Ymax1, Zmax1);
1882   B2.Get(Xmin2, Ymin2, Zmin2, Xmax2, Ymax2, Zmax2);
1883   if (Xmax2 < Xmin1 || Xmax1 < Xmin2 ||
1884       Ymax2 < Ymin1 || Ymax1 < Ymin2 ||
1885       Zmax2 < Zmin1 || Zmax1 < Zmin2) {
1886     return REL_NOT_CONNECTED;
1887   }
1888
1889   BRepExtrema_DistShapeShape dst (theBlock1, theBlock2);
1890   if (!dst.IsDone()) {
1891     return REL_UNKNOWN;
1892   }
1893
1894   if (dst.Value() > Precision::Confusion()) {
1895     return REL_NOT_CONNECTED;
1896   }
1897
1898   if (dst.InnerSolution()) {
1899     return REL_COLLISION_VV;
1900   }
1901
1902   Standard_Integer nbSol = dst.NbSolution();
1903   Standard_Integer relation = REL_OK;
1904   Standard_Integer nbVerts = 0;
1905   Standard_Integer nbEdges = 0;
1906   Standard_Integer sol = 1;
1907   for (; sol <= nbSol; sol++) {
1908     BRepExtrema_SupportType supp1 = dst.SupportTypeShape1(sol);
1909     BRepExtrema_SupportType supp2 = dst.SupportTypeShape2(sol);
1910     if (supp1 == BRepExtrema_IsVertex && supp2 == BRepExtrema_IsVertex) {
1911       nbVerts++;
1912     } else if (supp1 == BRepExtrema_IsInFace || supp2 == BRepExtrema_IsInFace) {
1913       return REL_COLLISION_FF;
1914     } else if (supp1 == BRepExtrema_IsOnEdge && supp2 == BRepExtrema_IsOnEdge) {
1915       nbEdges++;
1916     } else if ((supp1 == BRepExtrema_IsOnEdge && supp2 == BRepExtrema_IsVertex) ||
1917                (supp2 == BRepExtrema_IsOnEdge && supp1 == BRepExtrema_IsVertex)) {
1918       relation = REL_COLLISION_EE;
1919     } else {
1920     }
1921   }
1922
1923   if (relation != REL_OK) {
1924     return relation;
1925   }
1926
1927   TColStd_Array1OfInteger vertSol (1, nbVerts);
1928   TopTools_Array1OfShape V1 (1, nbVerts);
1929   TopTools_Array1OfShape V2 (1, nbVerts);
1930   Standard_Integer ivs = 0;
1931   for (sol = 1; sol <= nbSol; sol++) {
1932     if (dst.SupportTypeShape1(sol) == BRepExtrema_IsVertex &&
1933         dst.SupportTypeShape2(sol) == BRepExtrema_IsVertex) {
1934       TopoDS_Vertex Vcur = TopoDS::Vertex(dst.SupportOnShape1(sol));
1935       // Check, that this vertex is far enough from other solution vertices.
1936       Standard_Integer ii = 1;
1937       for (; ii <= ivs; ii++) {
1938         if (BRepTools::Compare(TopoDS::Vertex(V1(ii)), Vcur)) {
1939           continue;
1940         }
1941       }
1942       ivs++;
1943       vertSol(ivs) = sol;
1944       V1(ivs) = Vcur;
1945       V2(ivs) = dst.SupportOnShape2(sol);
1946     }
1947   }
1948
1949   // As we deal only with quadrangles,
1950   // 2, 3 or 4 vertex solutions can be found.
1951   if (ivs <= 1) {
1952     if (nbEdges > 0) {
1953       return REL_COLLISION_FF;
1954     }
1955     return REL_NOT_CONNECTED;
1956   }
1957   if (ivs > 4) {
1958     return REL_UNKNOWN;
1959   }
1960
1961   // Check sharing of coincident entities.
1962   if (ivs == 2 || ivs == 3) {
1963     // Map vertices and edges of the blocks
1964     TopTools_IndexedDataMapOfShapeListOfShape MVE1, MVE2;
1965     GEOMImpl_Block6Explorer::MapShapesAndAncestors
1966       (theBlock1, TopAbs_VERTEX, TopAbs_EDGE, MVE1);
1967     GEOMImpl_Block6Explorer::MapShapesAndAncestors
1968       (theBlock2, TopAbs_VERTEX, TopAbs_EDGE, MVE2);
1969
1970     if (ivs == 2) {
1971       // Find common edge
1972       TopoDS_Shape anEdge1, anEdge2;
1973       GEOMImpl_Block6Explorer::FindEdge(anEdge1, V1(1), V1(2), MVE1);
1974       if (anEdge1.IsNull()) return REL_UNKNOWN;
1975
1976       GEOMImpl_Block6Explorer::FindEdge(anEdge2, V2(1), V2(2), MVE2);
1977       if (anEdge2.IsNull()) return REL_UNKNOWN;
1978
1979       if (!anEdge1.IsSame(anEdge2)) return REL_NOT_GLUED;
1980
1981     } else { // ivs == 3
1982       // Find common edges
1983       Standard_Integer e1_v1 = 1;
1984       Standard_Integer e1_v2 = 2;
1985       Standard_Integer e2_v1 = 3;
1986       Standard_Integer e2_v2 = 1;
1987
1988       TopoDS_Shape anEdge11, anEdge12;
1989       GEOMImpl_Block6Explorer::FindEdge(anEdge11, V1(e1_v1), V1(e1_v2), MVE1);
1990       if (anEdge11.IsNull()) {
1991         e1_v2 = 3;
1992         e2_v1 = 2;
1993         GEOMImpl_Block6Explorer::FindEdge(anEdge11, V1(e1_v1), V1(e1_v2), MVE1);
1994         if (anEdge11.IsNull()) return REL_UNKNOWN;
1995       }
1996       GEOMImpl_Block6Explorer::FindEdge(anEdge12, V1(e2_v1), V1(e2_v2), MVE1);
1997       if (anEdge12.IsNull()) {
1998         e2_v2 = 5 - e2_v1;
1999         GEOMImpl_Block6Explorer::FindEdge(anEdge12, V1(e2_v1), V1(e2_v2), MVE1);
2000         if (anEdge12.IsNull()) return REL_UNKNOWN;
2001       }
2002
2003       TopoDS_Shape anEdge21, anEdge22;
2004       GEOMImpl_Block6Explorer::FindEdge(anEdge21, V2(e1_v1), V2(e1_v2), MVE2);
2005       if (anEdge21.IsNull()) return REL_UNKNOWN;
2006       GEOMImpl_Block6Explorer::FindEdge(anEdge22, V2(e2_v1), V2(e2_v2), MVE2);
2007       if (anEdge22.IsNull()) return REL_UNKNOWN;
2008
2009       // Check of edges coincidence (with some precision) have to be done here
2010       // if (!anEdge11.IsEqual(anEdge21)) return REL_UNKNOWN;
2011       // if (!anEdge12.IsEqual(anEdge22)) return REL_UNKNOWN;
2012
2013       // Check of edges sharing
2014       if (!anEdge11.IsSame(anEdge21)) return REL_NOT_GLUED;
2015       if (!anEdge12.IsSame(anEdge22)) return REL_NOT_GLUED;
2016     }
2017   }
2018
2019   if (ivs == 4) {
2020     // Map vertices and faces of the blocks
2021     TopTools_IndexedDataMapOfShapeListOfShape MVF1, MVF2;
2022     GEOMImpl_Block6Explorer::MapShapesAndAncestors
2023       (theBlock1, TopAbs_VERTEX, TopAbs_FACE, MVF1);
2024     GEOMImpl_Block6Explorer::MapShapesAndAncestors
2025       (theBlock2, TopAbs_VERTEX, TopAbs_FACE, MVF2);
2026
2027     TopoDS_Shape aFace1, aFace2;
2028     GEOMImpl_Block6Explorer::FindFace(aFace1, V1(1), V1(2), V1(3), V1(4), MVF1);
2029     if (aFace1.IsNull()) return REL_UNKNOWN;
2030     GEOMImpl_Block6Explorer::FindFace(aFace2, V2(1), V2(2), V2(3), V2(4), MVF2);
2031     if (aFace2.IsNull()) return REL_UNKNOWN;
2032
2033     // Check of faces coincidence (with some precision) have to be done here
2034     // if (!aFace1.IsEqual(aFace2)) return REL_UNKNOWN;
2035
2036     // Check of faces sharing
2037     if (!aFace1.IsSame(aFace2)) return REL_NOT_GLUED;
2038   }
2039
2040   return REL_OK;
2041 }
2042
2043 void FindConnected (const Standard_Integer         theBlockIndex,
2044                     const TColStd_Array2OfInteger& theRelations,
2045                     TColStd_MapOfInteger&          theProcessedMap,
2046                     TColStd_MapOfInteger&          theConnectedMap)
2047 {
2048   theConnectedMap.Add(theBlockIndex);
2049   theProcessedMap.Add(theBlockIndex);
2050
2051   Standard_Integer nbBlocks = theRelations.ColLength();
2052   Standard_Integer col = 1;
2053   for (; col <= nbBlocks; col++) {
2054     if (theRelations(theBlockIndex, col) == REL_OK ||
2055         theRelations(theBlockIndex, col) == REL_NOT_GLUED) {
2056       if (!theProcessedMap.Contains(col)) {
2057         FindConnected(col, theRelations, theProcessedMap, theConnectedMap);
2058       }
2059     }
2060   }
2061 }
2062
2063 Standard_Boolean HasAnyConnection (const Standard_Integer         theBlockIndex,
2064                                    const TColStd_MapOfInteger&    theWith,
2065                                    const TColStd_Array2OfInteger& theRelations,
2066                                    TColStd_MapOfInteger&          theProcessedMap)
2067 {
2068   theProcessedMap.Add(theBlockIndex);
2069
2070   Standard_Integer nbBlocks = theRelations.ColLength();
2071   Standard_Integer col = 1;
2072   for (; col <= nbBlocks; col++) {
2073     if (theRelations(theBlockIndex, col) != REL_NOT_CONNECTED) {
2074       if (!theProcessedMap.Contains(col)) {
2075         if (theWith.Contains(col))
2076           return Standard_True;
2077         if (HasAnyConnection(col, theWith, theRelations, theProcessedMap))
2078           return Standard_True;
2079       }
2080     }
2081   }
2082
2083   return Standard_False;
2084 }
2085
2086 //=============================================================================
2087 /*!
2088  *  PrintBCErrors
2089  */
2090 //=============================================================================
2091 TCollection_AsciiString GEOMImpl_IBlocksOperations::PrintBCErrors
2092                                                 (Handle(GEOM_Object)  theCompound,
2093                                                  const std::list<BCError>& theErrors)
2094 {
2095   TCollection_AsciiString aDescr;
2096
2097   std::list<BCError>::const_iterator errIt = theErrors.begin();
2098   int i = 0;
2099   for (; errIt != theErrors.end(); i++, errIt++) {
2100     BCError errStruct = *errIt;
2101
2102     switch (errStruct.error) {
2103     case NOT_BLOCK:
2104       aDescr += "\n\tNot a Blocks: ";
2105       break;
2106     case EXTRA_EDGE:
2107       aDescr += "\n\tHexahedral solids with degenerated and/or seam edges: ";
2108       break;
2109     case INVALID_CONNECTION:
2110       aDescr += "\n\tInvalid connection between two blocks: ";
2111       break;
2112     case NOT_CONNECTED:
2113       aDescr += "\n\tBlocks, not connected with main body: ";
2114       break;
2115     case NOT_GLUED:
2116       aDescr += "\n\tNot glued blocks: ";
2117       break;
2118     default:
2119       break;
2120     }
2121
2122     std::list<int> sshList = errStruct.incriminated;
2123     std::list<int>::iterator sshIt = sshList.begin();
2124     int jj = 0;
2125     for (; sshIt != sshList.end(); jj++, sshIt++) {
2126       if (jj > 0)
2127         aDescr += ", ";
2128       aDescr += TCollection_AsciiString(*sshIt);
2129     }
2130   }
2131
2132   return aDescr;
2133 }
2134
2135 //=============================================================================
2136 /*!
2137  *  CheckCompoundOfBlocks
2138  */
2139 //=============================================================================
2140 Standard_Boolean GEOMImpl_IBlocksOperations::CheckCompoundOfBlocks
2141                                               (Handle(GEOM_Object) theCompound,
2142                                                const Standard_Real theToleranceC1,
2143                                                std::list<BCError>& theErrors)
2144 {
2145   SetErrorCode(KO);
2146
2147   if (theCompound.IsNull()) return Standard_False;
2148   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2149
2150   Standard_Boolean isCompOfBlocks = Standard_True;
2151
2152   // Map sub-shapes and their indices
2153   TopTools_IndexedMapOfShape anIndices;
2154   TopExp::MapShapes(aBlockOrComp, anIndices);
2155
2156   // 1. Separate blocks from non-blocks
2157   TopTools_ListOfShape NOT; // Not blocks
2158   TopTools_ListOfShape EXT; // Hexahedral solids, having degenerated and/or seam edges
2159   TopTools_ListOfShape BLO; // All blocks from the given compound
2160   TopTools_ListOfShape NOQ; // All non-quadrangular faces
2161   AddBlocksFrom(aBlockOrComp, BLO, NOT, EXT, NOQ, theToleranceC1);
2162
2163   // Report non-blocks
2164   if (NOT.Extent() > 0) {
2165     isCompOfBlocks = Standard_False;
2166     BCError anErr;
2167     anErr.error = NOT_BLOCK;
2168     TopTools_ListIteratorOfListOfShape it (NOT);
2169     for (; it.More(); it.Next()) {
2170       anErr.incriminated.push_back(anIndices.FindIndex(it.Value()));
2171     }
2172     theErrors.push_back(anErr);
2173   }
2174
2175   // Report solids, having degenerated and/or seam edges
2176   if (EXT.Extent() > 0) {
2177     isCompOfBlocks = Standard_False;
2178     BCError anErr;
2179     anErr.error = EXTRA_EDGE;
2180     TopTools_ListIteratorOfListOfShape it (EXT);
2181     for (; it.More(); it.Next()) {
2182       anErr.incriminated.push_back(anIndices.FindIndex(it.Value()));
2183     }
2184     theErrors.push_back(anErr);
2185   }
2186
2187   Standard_Integer nbBlocks = BLO.Extent();
2188   if (nbBlocks == 0) {
2189     isCompOfBlocks = Standard_False;
2190     SetErrorCode(OK);
2191     return isCompOfBlocks;
2192   }
2193   if (nbBlocks == 1) {
2194     SetErrorCode(OK);
2195     return isCompOfBlocks;
2196   }
2197
2198   // Prepare data for 2. and 3.
2199   TColStd_Array2OfInteger aRelations (1, nbBlocks, 1, nbBlocks);
2200   aRelations.Init(REL_NOT_CONNECTED);
2201
2202   TopTools_IndexedMapOfShape mapBlocks;
2203
2204   BRep_Builder BB;
2205   TopoDS_Compound aComp;
2206   BB.MakeCompound(aComp);
2207
2208   TopTools_ListIteratorOfListOfShape BLOit (BLO);
2209   for (; BLOit.More(); BLOit.Next()) {
2210     mapBlocks.Add(BLOit.Value());
2211     BB.Add(aComp, BLOit.Value());
2212   }
2213
2214   // 2. Find glued blocks (having shared faces)
2215   TopTools_IndexedDataMapOfShapeListOfShape mapFaceBlocks;
2216   GEOMImpl_Block6Explorer::MapShapesAndAncestors
2217     (aComp, TopAbs_FACE, TopAbs_SOLID, mapFaceBlocks);
2218
2219   Standard_Integer prevInd = 0, curInd = 0;
2220   Standard_Integer ind = 1, nbFaces = mapFaceBlocks.Extent();
2221   for (; ind <= nbFaces; ind++) {
2222     const TopTools_ListOfShape& aGluedBlocks = mapFaceBlocks.FindFromIndex(ind);
2223     if (aGluedBlocks.Extent() > 1) { // Shared face found
2224       TopTools_ListIteratorOfListOfShape aGluedBlocksIt (aGluedBlocks);
2225       TopoDS_Shape prevBlock, curBlock;
2226       for (; aGluedBlocksIt.More(); aGluedBlocksIt.Next()) {
2227         curBlock = aGluedBlocksIt.Value();
2228         if (!prevBlock.IsNull()) {
2229           prevInd = mapBlocks.FindIndex(prevBlock);
2230           curInd  = mapBlocks.FindIndex(curBlock);
2231           aRelations.SetValue(prevInd, curInd, REL_OK);
2232           aRelations.SetValue(curInd, prevInd, REL_OK);
2233         }
2234         prevBlock = curBlock;
2235       }
2236     }
2237   }
2238
2239   // 3. Find not glued blocks
2240   GEOMAlgo_GlueAnalyser aGD;
2241
2242   aGD.SetShape(aComp);
2243   aGD.SetTolerance(Precision::Confusion());
2244   aGD.SetCheckGeometry(Standard_True);
2245   aGD.Perform();
2246
2247   Standard_Integer iErr, iWrn;
2248   iErr = aGD.ErrorStatus();
2249   if (iErr) {
2250     SetErrorCode("Error in GEOMAlgo_GlueAnalyser");
2251     return isCompOfBlocks;
2252   }
2253   iWrn = aGD.WarningStatus();
2254   if (iWrn) {
2255     MESSAGE("Warning in GEOMAlgo_GlueAnalyser");
2256   }
2257
2258   // Report not glued blocks
2259   if (aGD.HasSolidsToGlue()) {
2260     isCompOfBlocks = Standard_False;
2261     Standard_Integer aSx1Ind, aSx2Ind;
2262
2263     const GEOMAlgo_ListOfCoupleOfShapes& aLCS = aGD.SolidsToGlue();
2264     GEOMAlgo_ListIteratorOfListOfCoupleOfShapes aItCS (aLCS);
2265     for (; aItCS.More(); aItCS.Next()) {
2266       const GEOMAlgo_CoupleOfShapes& aCS = aItCS.Value();
2267       const TopoDS_Shape& aSx1 = aCS.Shape1();
2268       const TopoDS_Shape& aSx2 = aCS.Shape2();
2269
2270       aSx1Ind = mapBlocks.FindIndex(aSx1);
2271       aSx2Ind = mapBlocks.FindIndex(aSx2);
2272       aRelations.SetValue(aSx1Ind, aSx2Ind, NOT_GLUED);
2273       aRelations.SetValue(aSx2Ind, aSx1Ind, NOT_GLUED);
2274
2275       BCError anErr;
2276       anErr.error = NOT_GLUED;
2277       anErr.incriminated.push_back(anIndices.FindIndex(aSx1));
2278       anErr.incriminated.push_back(anIndices.FindIndex(aSx2));
2279       theErrors.push_back(anErr);
2280     }
2281   }
2282
2283   // 4. Find largest set of connected (good connection or not glued) blocks
2284   Standard_Integer ibl = 1;
2285   TColStd_MapOfInteger aProcessedMap;
2286   TColStd_MapOfInteger aLargestSet;
2287   TColStd_MapOfInteger aCurrentSet;
2288   for (ibl = 1; ibl <= nbBlocks; ibl++) {
2289     if (!aProcessedMap.Contains(ibl)) {
2290       aCurrentSet.Clear();
2291       FindConnected(ibl, aRelations, aProcessedMap, aCurrentSet);
2292       if (aCurrentSet.Extent() > aLargestSet.Extent()) {
2293         aLargestSet = aCurrentSet;
2294       }
2295     }
2296   }
2297
2298   // 5. Report all blocks, isolated from <aLargestSet>
2299   BCError anErr;
2300   anErr.error = NOT_CONNECTED;
2301   Standard_Boolean hasIsolated = Standard_False;
2302   for (ibl = 1; ibl <= nbBlocks; ibl++) {
2303     if (!aLargestSet.Contains(ibl)) {
2304       aProcessedMap.Clear();
2305       if (!HasAnyConnection(ibl, aLargestSet, aRelations, aProcessedMap)) {
2306         // report connection absence
2307         hasIsolated = Standard_True;
2308         anErr.incriminated.push_back(anIndices.FindIndex(mapBlocks.FindKey(ibl)));
2309       }
2310     }
2311   }
2312   if (hasIsolated) {
2313     isCompOfBlocks = Standard_False;
2314     theErrors.push_back(anErr);
2315   }
2316
2317   SetErrorCode(OK);
2318   return isCompOfBlocks;
2319 }
2320
2321 //=============================================================================
2322 /*!
2323  *  GetNonBlocks
2324  */
2325 //=============================================================================
2326 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetNonBlocks
2327                                      (Handle(GEOM_Object)  theShape,
2328                                       const Standard_Real  theToleranceC1,
2329                                       Handle(GEOM_Object)& theNonQuads)
2330 {
2331   SetErrorCode(KO);
2332
2333   if (theShape.IsNull()) return NULL;
2334   TopoDS_Shape aShape = theShape->GetValue();
2335
2336   // Separate blocks from non-blocks
2337   TopTools_ListOfShape BLO; // All blocks from the given compound
2338   TopTools_ListOfShape NOT; // Not blocks
2339   TopTools_ListOfShape EXT; // Hexahedral solids, having degenerated and/or seam edges
2340   TopTools_ListOfShape NOQ; // All non-quadrangular faces
2341   AddBlocksFrom(aShape, BLO, NOT, EXT, NOQ, theToleranceC1);
2342
2343   if (NOT.IsEmpty() && EXT.IsEmpty() && NOQ.IsEmpty()) {
2344     SetErrorCode("NOT_FOUND_ANY");
2345     return NULL;
2346   }
2347
2348   // Map sub-shapes and their indices
2349   TopTools_IndexedMapOfShape anIndices;
2350   TopExp::MapShapes(aShape, anIndices);
2351
2352   // Non-blocks
2353   Handle(GEOM_Object) aNonBlocks;
2354   if (NOT.Extent() > 0 || EXT.Extent() > 0) {
2355     Handle(TColStd_HArray1OfInteger) anArray =
2356       new TColStd_HArray1OfInteger (1, NOT.Extent() + EXT.Extent());
2357     Standard_Integer ii = 1;
2358     TopTools_ListIteratorOfListOfShape it1 (NOT);
2359     for (; it1.More(); it1.Next(), ii++) {
2360       anArray->SetValue(ii, anIndices.FindIndex(it1.Value()));
2361     }
2362     TopTools_ListIteratorOfListOfShape it2 (EXT);
2363     for (; it2.More(); it2.Next(), ii++) {
2364       anArray->SetValue(ii, anIndices.FindIndex(it2.Value()));
2365     }
2366
2367     aNonBlocks = GetEngine()->AddSubShape(theShape, anArray);
2368     if (aNonBlocks.IsNull()) {
2369       SetErrorCode("Error in algorithm: result found, but cannot be returned.");
2370       return NULL;
2371     }
2372     aNonBlocks->SetType(GEOM_GROUP);
2373     TDF_Label aFreeLabel = aNonBlocks->GetFreeLabel();
2374     TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)TopAbs_SOLID);
2375   }
2376
2377   // Non-quadrangles
2378   if (NOQ.Extent() > 0) {
2379     Handle(TColStd_HArray1OfInteger) anArray =
2380       new TColStd_HArray1OfInteger (1, NOQ.Extent());
2381     Standard_Integer ii = 1;
2382     TopTools_ListIteratorOfListOfShape it1 (NOQ);
2383     for (; it1.More(); it1.Next(), ii++) {
2384       anArray->SetValue(ii, anIndices.FindIndex(it1.Value()));
2385     }
2386
2387     theNonQuads = GetEngine()->AddSubShape(theShape, anArray);
2388     if (theNonQuads.IsNull()) {
2389       SetErrorCode("Error in algorithm: result found, but cannot be returned.");
2390       return NULL;
2391     }
2392     theNonQuads->SetType(GEOM_GROUP);
2393     TDF_Label aFreeLabel = theNonQuads->GetFreeLabel();
2394     TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)TopAbs_FACE);
2395   }
2396
2397   //Make a Python command
2398   Handle(GEOM_Function) aMainShape = theShape->GetLastFunction();
2399   GEOM::TPythonDump pd (aMainShape, /*append=*/true);
2400   pd << "(";
2401   if (aNonBlocks.IsNull())
2402     pd << "no_bad_solids";
2403   else
2404     pd << aNonBlocks;
2405   pd << ", ";
2406   if (theNonQuads.IsNull())
2407     pd << "no_bad_faces";
2408   else
2409     pd << theNonQuads;
2410   pd << ") = geompy.GetNonBlocks(" << theShape << ")";
2411
2412   SetErrorCode(OK);
2413   return aNonBlocks;
2414 }
2415
2416 //=============================================================================
2417 /*!
2418  *  RemoveExtraEdges
2419  */
2420 //=============================================================================
2421 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::RemoveExtraEdges
2422                                      (Handle(GEOM_Object) theObject,
2423                                       const Standard_Integer theOptimumNbFaces)
2424 {
2425   SetErrorCode(KO);
2426
2427   if (theObject.IsNull()) return NULL;
2428
2429   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
2430   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be fixed
2431
2432   //Add a new Copy object
2433   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
2434
2435   //Add a function
2436   Handle(GEOM_Function) aFunction =
2437     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_REMOVE_EXTRA);
2438
2439   //Check if the function is set correctly
2440   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
2441
2442   GEOMImpl_IBlockTrsf aTI (aFunction);
2443   aTI.SetOriginal(aLastFunction);
2444   aTI.SetOptimumNbFaces(theOptimumNbFaces);
2445
2446   //Compute the fixed shape
2447   try {
2448     OCC_CATCH_SIGNALS;
2449     if (!GetSolver()->ComputeFunction(aFunction)) {
2450       SetErrorCode("Block driver failed to remove extra edges of the given shape");
2451       return NULL;
2452     }
2453   }
2454   catch (Standard_Failure) {
2455     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2456     SetErrorCode(aFail->GetMessageString());
2457     return NULL;
2458   }
2459
2460   //Make a Python command
2461   std::string doUnionFaces = (theOptimumNbFaces < 0) ? "False" : "True";
2462   GEOM::TPythonDump(aFunction) << aCopy << " = geompy.RemoveExtraEdges("
2463                                << theObject << ", " << doUnionFaces.data() << ")";
2464
2465   SetErrorCode(OK);
2466   return aCopy;
2467 }
2468
2469 //=============================================================================
2470 /*!
2471  *  UnionFaces
2472  */
2473 //=============================================================================
2474 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::UnionFaces
2475                                      (Handle(GEOM_Object) theObject)
2476 {
2477   SetErrorCode(KO);
2478
2479   if (theObject.IsNull()) return NULL;
2480
2481   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
2482   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be fixed
2483
2484   //Add a new Copy object
2485   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
2486
2487   //Add a function
2488   Handle(GEOM_Function) aFunction =
2489     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_UNION_FACES);
2490
2491   //Check if the function is set correctly
2492   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
2493
2494   GEOMImpl_IBlockTrsf aTI (aFunction);
2495   aTI.SetOriginal(aLastFunction);
2496
2497   //Compute the fixed shape
2498   try {
2499     OCC_CATCH_SIGNALS;
2500     if (!GetSolver()->ComputeFunction(aFunction)) {
2501       SetErrorCode("Block driver failed to remove extra edges of the given shape");
2502       return NULL;
2503     }
2504   }
2505   catch (Standard_Failure) {
2506     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2507     SetErrorCode(aFail->GetMessageString());
2508     return NULL;
2509   }
2510
2511   //Make a Python command
2512   GEOM::TPythonDump(aFunction) << aCopy << " = geompy.UnionFaces("
2513                                << theObject << ")";
2514
2515   SetErrorCode(OK);
2516   return aCopy;
2517 }
2518
2519 //=============================================================================
2520 /*!
2521  *  CheckAndImprove
2522  */
2523 //=============================================================================
2524 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::CheckAndImprove
2525                                              (Handle(GEOM_Object) theObject)
2526 {
2527   SetErrorCode(KO);
2528
2529   if (theObject.IsNull()) return NULL;
2530
2531   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
2532   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be fixed
2533
2534   //Add a new Copy object
2535   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
2536
2537   //Add a function
2538   Handle(GEOM_Function) aFunction =
2539     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_COMPOUND_IMPROVE);
2540
2541   //Check if the function is set correctly
2542   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
2543
2544   GEOMImpl_IBlockTrsf aTI (aFunction);
2545   aTI.SetOriginal(aLastFunction);
2546
2547   // -1 means do not unite faces on common surface (?except case of seam edge between them?)
2548   //aTI.SetOptimumNbFaces(-1);
2549   aTI.SetOptimumNbFaces(6);
2550
2551   //Compute the fixed shape
2552   try {
2553     OCC_CATCH_SIGNALS;
2554     if (!GetSolver()->ComputeFunction(aFunction)) {
2555       SetErrorCode("Block driver failed to improve the given blocks compound");
2556       return NULL;
2557     }
2558   }
2559   catch (Standard_Failure) {
2560     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2561     SetErrorCode(aFail->GetMessageString());
2562     return NULL;
2563   }
2564
2565   //Make a Python command
2566   GEOM::TPythonDump(aFunction) << aCopy
2567     << " = geompy.CheckAndImprove(" << theObject << ")";
2568
2569   SetErrorCode(OK);
2570   return aCopy;
2571 }
2572
2573 //=============================================================================
2574 /*!
2575  *  ExplodeCompoundOfBlocks
2576  */
2577 //=============================================================================
2578 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::ExplodeCompoundOfBlocks
2579                                                 (Handle(GEOM_Object)    theCompound,
2580                                                  const Standard_Integer theMinNbFaces,
2581                                                  const Standard_Integer theMaxNbFaces)
2582 {
2583   SetErrorCode(KO);
2584
2585   if (theCompound.IsNull()) return NULL;
2586   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2587   if (aBlockOrComp.IsNull()) return NULL;
2588
2589   Handle(TColStd_HSequenceOfTransient) aBlocks = new TColStd_HSequenceOfTransient;
2590   Handle(GEOM_Object) anObj;
2591   Handle(GEOM_Function) aFunction;
2592
2593   TopTools_MapOfShape mapShape;
2594   TCollection_AsciiString anAsciiList, anEntry;
2595
2596   // Map shapes
2597   TopTools_IndexedMapOfShape anIndices;
2598   TopExp::MapShapes(aBlockOrComp, anIndices);
2599   Handle(TColStd_HArray1OfInteger) anArray;
2600
2601   // Explode
2602   try {
2603     OCC_CATCH_SIGNALS;
2604     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
2605     for (; exp.More(); exp.Next()) {
2606       if (mapShape.Add(exp.Current())) {
2607         TopoDS_Shape aSolid = exp.Current();
2608
2609         TopTools_MapOfShape mapFaces;
2610         TopExp_Explorer expF (aSolid, TopAbs_FACE);
2611         Standard_Integer nbFaces = 0;
2612         for (; expF.More(); expF.Next()) {
2613           if (mapFaces.Add(expF.Current())) {
2614             nbFaces++;
2615           }
2616         }
2617
2618         if (theMinNbFaces <= nbFaces && nbFaces <= theMaxNbFaces) {
2619           anArray = new TColStd_HArray1OfInteger(1,1);
2620           anArray->SetValue(1, anIndices.FindIndex(aSolid));
2621           anObj = GetEngine()->AddSubShape(theCompound, anArray);
2622           aBlocks->Append(anObj);
2623
2624           //Make a Python command
2625           TDF_Tool::Entry(anObj->GetEntry(), anEntry);
2626           anAsciiList += anEntry + ", ";
2627         }
2628       }
2629     }
2630   }
2631   catch (Standard_Failure) {
2632     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2633     SetErrorCode(aFail->GetMessageString());
2634     return aBlocks;
2635   }
2636
2637   if (aBlocks->IsEmpty()) {
2638     SetErrorCode("There are no specified blocks in the given shape");
2639     return aBlocks;
2640   }
2641
2642   anAsciiList.Trunc(anAsciiList.Length() - 2);
2643
2644   //The explode doesn't change object so no new function is required.
2645   aFunction = theCompound->GetLastFunction();
2646
2647   //Make a Python command
2648   GEOM::TPythonDump(aFunction, /*append=*/true)
2649     << "[" << anAsciiList.ToCString() << "] = geompy.MakeBlockExplode("
2650     << theCompound << ", " << theMinNbFaces << ", " << theMaxNbFaces << ")";
2651
2652   SetErrorCode(OK);
2653   return aBlocks;
2654 }
2655
2656 //=============================================================================
2657 /*!
2658  *  GetBlockNearPoint
2659  */
2660 //=============================================================================
2661 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockNearPoint
2662                                                 (Handle(GEOM_Object) theCompound,
2663                                                  Handle(GEOM_Object) thePoint)
2664 {
2665   SetErrorCode(KO);
2666
2667   //New object
2668   Handle(GEOM_Object) aResult;
2669
2670   // Arguments
2671   if (theCompound.IsNull() || thePoint.IsNull()) return NULL;
2672
2673   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2674   if (aBlockOrComp.IsNull()) {
2675     SetErrorCode("Compound is null");
2676     return NULL;
2677   }
2678   if (aBlockOrComp.ShapeType() != TopAbs_COMPOUND &&
2679       aBlockOrComp.ShapeType() != TopAbs_COMPSOLID) {
2680     SetErrorCode("Shape to find block in is not a compound");
2681     return NULL;
2682   }
2683
2684   TopoDS_Shape anArg = thePoint->GetValue();
2685   if (anArg.IsNull()) {
2686     SetErrorCode("Point is null");
2687     return NULL;
2688   }
2689   if (anArg.ShapeType() != TopAbs_VERTEX) {
2690     SetErrorCode("Shape for block identification is not a vertex");
2691     return NULL;
2692   }
2693
2694   //Compute the Block value
2695   try {
2696     OCC_CATCH_SIGNALS;
2697     TopoDS_Shape aShape;
2698
2699     TopoDS_Vertex aVert = TopoDS::Vertex(anArg);
2700     gp_Pnt aPnt = BRep_Tool::Pnt(aVert);
2701     Standard_Real PX, PY, PZ;
2702     aPnt.Coord(PX, PY, PZ);
2703
2704     // 1. Classify the point relatively each block
2705     Standard_Integer nearest = 2, nbFound = 0;
2706     TopTools_DataMapOfShapeInteger mapShapeDist;
2707     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
2708     for (; exp.More(); exp.Next()) {
2709       TopoDS_Shape aSolid = exp.Current();
2710
2711       if (!mapShapeDist.IsBound(aSolid)) {
2712         Standard_Integer aDistance = 2;
2713
2714         // 1.a. Classify relatively Bounding box
2715         Standard_Real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax;
2716         Bnd_Box BB;
2717         BRepBndLib::Add(aSolid, BB);
2718         BB.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
2719         if (PX < Xmin || Xmax < PX ||
2720             PY < Ymin || Ymax < PY ||
2721             PZ < Zmin || Zmax < PZ) {
2722           // OUT of bounding box
2723           aDistance = 1;
2724         } else {
2725           // 1.b. Classify relatively the solid itself
2726           BRepClass3d_SolidClassifier SC (aSolid, aPnt, Precision::Confusion());
2727           if (SC.State() == TopAbs_IN) {
2728             aDistance = -1;
2729           } else if (SC.State() == TopAbs_ON) {
2730             aDistance = 0;
2731           } else { // OUT
2732             aDistance = 1;
2733           }
2734         }
2735
2736         if (aDistance < nearest) {
2737           nearest = aDistance;
2738           aShape = aSolid;
2739           nbFound = 1;
2740
2741           // A first found block, containing the point inside, will be returned.
2742           // It is the solution, if there are no intersecting blocks in the compound.
2743           if (nearest == -1) break;
2744
2745         } else if (aDistance == nearest) {
2746           nbFound++;
2747         } else {
2748         }
2749
2750         mapShapeDist.Bind(aSolid, aDistance);
2751       } // if (!mapShapeDist.IsBound(aSolid))
2752     }
2753
2754     // 2. Define block, containing the point or having minimum distance to it
2755     if (nbFound > 1) {
2756       if (nearest == 0) {
2757         // The point is on boundary of some blocks and there are
2758         // no blocks, having the point inside their volume
2759         SetErrorCode("Multiple blocks near the given point are found");
2760         return NULL;
2761
2762       } else if (nearest == 1) {
2763         // The point is outside some blocks and there are
2764         // no blocks, having the point inside or on boundary.
2765         // We will get a nearest block
2766         Standard_Real minDist = RealLast();
2767         TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
2768         for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
2769           if (mapShapeDistIter.Value() == 1) {
2770             TopoDS_Shape aSolid = mapShapeDistIter.Key();
2771             BRepExtrema_DistShapeShape aDistTool (aVert, aSolid);
2772             if (!aDistTool.IsDone()) {
2773               SetErrorCode("Can not find a distance from the given point to one of blocks");
2774               return NULL;
2775             }
2776             Standard_Real aDist = aDistTool.Value();
2777             if (aDist < minDist) {
2778               minDist = aDist;
2779               aShape = aSolid;
2780             }
2781           }
2782         }
2783       } else { // nearest == -1
2784 //        // The point is inside some blocks.
2785 //        // We will get a block with nearest center
2786 //        Standard_Real minDist = RealLast();
2787 //        TopTools_DataMapIteratorOfDataMapOfShapeInteger mapShapeDistIter (mapShapeDist);
2788 //        for (; mapShapeDistIter.More(); mapShapeDistIter.Next()) {
2789 //          if (mapShapeDistIter.Value() == -1) {
2790 //            TopoDS_Shape aSolid = mapShapeDistIter.Key();
2791 //            GProp_GProps aSystem;
2792 //            BRepGProp::VolumeProperties(aSolid, aSystem);
2793 //            gp_Pnt aCenterMass = aSystem.CentreOfMass();
2794 //
2795 //            Standard_Real aDist = aCenterMass.Distance(aPnt);
2796 //            if (aDist < minDist) {
2797 //              minDist = aDist;
2798 //              aShape = aSolid;
2799 //            }
2800 //          }
2801 //        }
2802       }
2803     } // if (nbFound > 1)
2804
2805     if (nbFound == 0) {
2806       SetErrorCode("There are no blocks near the given point");
2807       return NULL;
2808     } else {
2809       TopTools_IndexedMapOfShape anIndices;
2810       TopExp::MapShapes(aBlockOrComp, anIndices);
2811       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
2812       anArray->SetValue(1, anIndices.FindIndex(aShape));
2813       aResult = GetEngine()->AddSubShape(theCompound, anArray);
2814     }
2815   }
2816   catch (Standard_Failure) {
2817     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2818     SetErrorCode(aFail->GetMessageString());
2819     return NULL;
2820   }
2821
2822   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
2823
2824   //Make a Python command
2825   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetBlockNearPoint("
2826                                << theCompound << ", " << thePoint << ")";
2827
2828   SetErrorCode(OK);
2829   return aResult;
2830 }
2831
2832 //=============================================================================
2833 /*!
2834  *  GetBlockByParts
2835  */
2836 //=============================================================================
2837 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::GetBlockByParts
2838                       (Handle(GEOM_Object)                         theCompound,
2839                        const Handle(TColStd_HSequenceOfTransient)& theParts)
2840 {
2841   SetErrorCode(KO);
2842
2843   Handle(GEOM_Object) aResult;
2844
2845   if (theCompound.IsNull() || theParts.IsNull()) return NULL;
2846   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2847   if (aBlockOrComp.IsNull()) return NULL;
2848
2849   //Get the parts
2850   Standard_Integer argi, aLen = theParts->Length();
2851   TopTools_Array1OfShape anArgs (1, aLen);
2852   TCollection_AsciiString anEntry, aPartsDescr;
2853   for (argi = 1; argi <= aLen; argi++) {
2854     Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast(theParts->Value(argi));
2855     Handle(GEOM_Function) aRef = anObj->GetLastFunction();
2856     if (aRef.IsNull()) return NULL;
2857
2858     TopoDS_Shape anArg = aRef->GetValue();
2859     if (anArg.IsNull()) {
2860       SetErrorCode("Null shape is given as argument");
2861       return NULL;
2862     }
2863     anArgs(argi) = anArg;
2864
2865     // For Python command
2866     TDF_Tool::Entry(anObj->GetEntry(), anEntry);
2867     if (argi > 1) aPartsDescr += ", ";
2868     aPartsDescr += anEntry;
2869   }
2870
2871   //Compute the Block value
2872   try {
2873     OCC_CATCH_SIGNALS;
2874     // 1. Explode compound on solids
2875     TopTools_MapOfShape mapShape;
2876     Standard_Integer nbSolids = 0;
2877     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
2878     for (; exp.More(); exp.Next()) {
2879       if (mapShape.Add(exp.Current())) {
2880         nbSolids++;
2881       }
2882     }
2883
2884     mapShape.Clear();
2885     Standard_Integer ind = 1;
2886     TopTools_Array1OfShape aSolids (1, nbSolids);
2887     TColStd_Array1OfInteger aNbParts (1, nbSolids);
2888     for (exp.Init(aBlockOrComp, TopAbs_SOLID); exp.More(); exp.Next(), ind++) {
2889       if (mapShape.Add(exp.Current())) {
2890         TopoDS_Shape aSolid = exp.Current();
2891         aSolids(ind) = aSolid;
2892         aNbParts(ind) = 0;
2893
2894         // 2. Define quantity of parts, contained in each solid
2895         TopTools_IndexedMapOfShape aSubShapes;
2896         TopExp::MapShapes(aSolid, aSubShapes);
2897         for (argi = 1; argi <= aLen; argi++) {
2898           if (aSubShapes.Contains(anArgs(argi))) {
2899             aNbParts(ind)++;
2900           }
2901         }
2902       }
2903     }
2904
2905     // 3. Define solid, containing maximum quantity of parts
2906     Standard_Integer maxNb = 0, nbFound = 0;
2907     TopoDS_Shape aShape;
2908     for (ind = 1; ind <= nbSolids; ind++) {
2909       if (aNbParts(ind) > maxNb) {
2910         maxNb = aNbParts(ind);
2911         aShape = aSolids(ind);
2912         nbFound = 1;
2913       } else if (aNbParts(ind) == maxNb) {
2914         nbFound++;
2915       } else {
2916       }
2917     }
2918     if (nbFound > 1) {
2919       SetErrorCode("Multiple blocks, containing maximum quantity of the given parts, are found");
2920       return NULL;
2921     } else if (nbFound == 0) {
2922       SetErrorCode("There are no blocks, containing the given parts");
2923       return NULL;
2924     } else {
2925       TopTools_IndexedMapOfShape anIndices;
2926       TopExp::MapShapes(aBlockOrComp, anIndices);
2927       Handle(TColStd_HArray1OfInteger) anArray = new TColStd_HArray1OfInteger(1,1);
2928       anArray->SetValue(1, anIndices.FindIndex(aShape));
2929       aResult = GetEngine()->AddSubShape(theCompound, anArray);
2930     }
2931   } catch (Standard_Failure) {
2932     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2933     SetErrorCode(aFail->GetMessageString());
2934     return NULL;
2935   }
2936
2937   Handle(GEOM_Function) aFunction = aResult->GetLastFunction();
2938
2939   //Make a Python command
2940   GEOM::TPythonDump(aFunction) << aResult << " = geompy.GetBlockByParts("
2941     << theCompound << ", [" << aPartsDescr.ToCString() << "])";
2942
2943   SetErrorCode(OK);
2944   return aResult;
2945 }
2946
2947 //=============================================================================
2948 /*!
2949  *  GetBlocksByParts
2950  */
2951 //=============================================================================
2952 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::GetBlocksByParts
2953                       (Handle(GEOM_Object)                         theCompound,
2954                        const Handle(TColStd_HSequenceOfTransient)& theParts)
2955 {
2956   SetErrorCode(KO);
2957
2958   if (theCompound.IsNull() || theParts.IsNull()) return NULL;
2959   TopoDS_Shape aBlockOrComp = theCompound->GetValue();
2960   if (aBlockOrComp.IsNull()) return NULL;
2961
2962   Handle(TColStd_HSequenceOfTransient) aBlocks = new TColStd_HSequenceOfTransient;
2963   Handle(GEOM_Object) anObj;
2964   Handle(GEOM_Function) aFunction;
2965
2966   //Get the parts
2967   Standard_Integer argi, aLen = theParts->Length();
2968   TopTools_Array1OfShape anArgs (1, aLen);
2969   TCollection_AsciiString anEntry, aPartsDescr, anAsciiList;
2970
2971   for (argi = 1; argi <= aLen; argi++) {
2972     Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast(theParts->Value(argi));
2973     Handle(GEOM_Function) aRef = anObj->GetLastFunction();
2974     if (aRef.IsNull()) return NULL;
2975
2976     TopoDS_Shape anArg = aRef->GetValue();
2977     if (anArg.IsNull()) {
2978       SetErrorCode("Null shape is given as argument");
2979       return NULL;
2980     }
2981     anArgs(argi) = anArg;
2982
2983     // For Python command
2984     TDF_Tool::Entry(anObj->GetEntry(), anEntry);
2985     aPartsDescr += anEntry + ", ";
2986   }
2987
2988   //Get the Blocks
2989   try {
2990     OCC_CATCH_SIGNALS;
2991     TopTools_MapOfShape mapShape;
2992     Standard_Integer nbSolids = 0;
2993     TopExp_Explorer exp (aBlockOrComp, TopAbs_SOLID);
2994     for (; exp.More(); exp.Next()) {
2995       if (mapShape.Add(exp.Current())) {
2996         nbSolids++;
2997       }
2998     }
2999
3000     mapShape.Clear();
3001     Standard_Integer ind = 1;
3002     TopTools_Array1OfShape aSolids (1, nbSolids);
3003     TColStd_Array1OfInteger aNbParts (1, nbSolids);
3004     for (exp.Init(aBlockOrComp, TopAbs_SOLID); exp.More(); exp.Next(), ind++) {
3005       if (mapShape.Add(exp.Current())) {
3006         TopoDS_Shape aSolid = exp.Current();
3007         aSolids(ind) = aSolid;
3008         aNbParts(ind) = 0;
3009
3010         // 2. Define quantity of parts, contained in each solid
3011         TopTools_IndexedMapOfShape aSubShapes;
3012         TopExp::MapShapes(aSolid, aSubShapes);
3013         for (argi = 1; argi <= aLen; argi++) {
3014           if (aSubShapes.Contains(anArgs(argi))) {
3015             aNbParts(ind)++;
3016           }
3017         }
3018       }
3019     }
3020
3021     // 3. Define solid, containing maximum quantity of parts
3022     Standard_Integer maxNb = 0, nbFound = 0;
3023     for (ind = 1; ind <= nbSolids; ind++) {
3024       if (aNbParts(ind) > maxNb) {
3025         maxNb = aNbParts(ind);
3026         nbFound = 1;
3027       } else if (aNbParts(ind) == maxNb) {
3028         nbFound++;
3029       } else {
3030       }
3031     }
3032     if (nbFound == 0) {
3033       SetErrorCode("There are no blocks, containing the given parts");
3034       return NULL;
3035     }
3036
3037     // Map shapes
3038     TopTools_IndexedMapOfShape anIndices;
3039     TopExp::MapShapes(aBlockOrComp, anIndices);
3040     Handle(TColStd_HArray1OfInteger) anArray;
3041
3042     for (ind = 1; ind <= nbSolids; ind++) {
3043       if (aNbParts(ind) == maxNb) {
3044         anArray = new TColStd_HArray1OfInteger(1,1);
3045         anArray->SetValue(1, anIndices.FindIndex(aSolids(ind)));
3046         anObj = GetEngine()->AddSubShape(theCompound, anArray);
3047         aBlocks->Append(anObj);
3048
3049         // For Python command
3050         TDF_Tool::Entry(anObj->GetEntry(), anEntry);
3051         anAsciiList += anEntry + ", ";
3052         if (aFunction.IsNull())
3053           aFunction = anObj->GetLastFunction();
3054       }
3055     }
3056   }
3057   catch (Standard_Failure) {
3058     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3059     SetErrorCode(aFail->GetMessageString());
3060     return NULL;
3061   }
3062
3063   //Make a Python command
3064   aPartsDescr.Trunc(aPartsDescr.Length() - 2);
3065   anAsciiList.Trunc(anAsciiList.Length() - 2);
3066
3067   GEOM::TPythonDump(aFunction) << "[" << anAsciiList.ToCString()
3068     << "] = geompy.GetBlocksByParts(" << theCompound
3069       << ", [" << aPartsDescr.ToCString() << "])";
3070
3071   SetErrorCode(OK);
3072   return aBlocks;
3073 }
3074
3075 //=============================================================================
3076 /*!
3077  *  MakeMultiTransformation1D
3078  */
3079 //=============================================================================
3080 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation1D
3081                                              (Handle(GEOM_Object)    theObject,
3082                                               const Standard_Integer theDirFace1,
3083                                               const Standard_Integer theDirFace2,
3084                                               const Standard_Integer theNbTimes)
3085 {
3086   SetErrorCode(KO);
3087
3088   if (theObject.IsNull()) return NULL;
3089
3090   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
3091   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be moved
3092
3093   //Add a new Copy object
3094   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
3095
3096   //Add a translate function
3097   Handle(GEOM_Function) aFunction =
3098     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_MULTI_TRANSFORM_1D);
3099
3100   //Check if the function is set correctly
3101   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
3102
3103   GEOMImpl_IBlockTrsf aTI (aFunction);
3104   aTI.SetOriginal(aLastFunction);
3105   aTI.SetFace1U(theDirFace1);
3106   aTI.SetFace2U(theDirFace2);
3107   aTI.SetNbIterU(theNbTimes);
3108
3109   //Compute the transformation
3110   try {
3111     OCC_CATCH_SIGNALS;
3112     if (!GetSolver()->ComputeFunction(aFunction)) {
3113       SetErrorCode("Block driver failed to make multi-transformation");
3114       return NULL;
3115     }
3116   }
3117   catch (Standard_Failure) {
3118     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3119     SetErrorCode(aFail->GetMessageString());
3120     return NULL;
3121   }
3122
3123   //Make a Python command
3124   GEOM::TPythonDump(aFunction) << aCopy << " = geompy.MakeMultiTransformation1D("
3125     << theObject << ", " << theDirFace1 << ", " << theDirFace2 << ", " << theNbTimes << ")";
3126
3127   SetErrorCode(OK);
3128   return aCopy;
3129 }
3130
3131 //=============================================================================
3132 /*!
3133  *  MakeMultiTransformation2D
3134  */
3135 //=============================================================================
3136 Handle(GEOM_Object) GEOMImpl_IBlocksOperations::MakeMultiTransformation2D
3137                                              (Handle(GEOM_Object)    theObject,
3138                                               const Standard_Integer theDirFace1U,
3139                                               const Standard_Integer theDirFace2U,
3140                                               const Standard_Integer theNbTimesU,
3141                                               const Standard_Integer theDirFace1V,
3142                                               const Standard_Integer theDirFace2V,
3143                                               const Standard_Integer theNbTimesV)
3144 {
3145   SetErrorCode(KO);
3146
3147   if (theObject.IsNull()) return NULL;
3148
3149   Handle(GEOM_Function) aLastFunction = theObject->GetLastFunction();
3150   if (aLastFunction.IsNull()) return NULL; //There is no function which creates an object to be moved
3151
3152   //Add a new Copy object
3153   Handle(GEOM_Object) aCopy = GetEngine()->AddObject(GetDocID(), GEOM_COPY);
3154
3155   //Add a translate function
3156   Handle(GEOM_Function) aFunction =
3157     aCopy->AddFunction(GEOMImpl_BlockDriver::GetID(), BLOCK_MULTI_TRANSFORM_2D);
3158
3159   //Check if the function is set correctly
3160   if (aFunction->GetDriverGUID() != GEOMImpl_BlockDriver::GetID()) return NULL;
3161
3162   GEOMImpl_IBlockTrsf aTI (aFunction);
3163   aTI.SetOriginal(aLastFunction);
3164   aTI.SetFace1U(theDirFace1U);
3165   aTI.SetFace2U(theDirFace2U);
3166   aTI.SetNbIterU(theNbTimesU);
3167   aTI.SetFace1V(theDirFace1V);
3168   aTI.SetFace2V(theDirFace2V);
3169   aTI.SetNbIterV(theNbTimesV);
3170
3171   //Compute the transformation
3172   try {
3173     OCC_CATCH_SIGNALS;
3174     if (!GetSolver()->ComputeFunction(aFunction)) {
3175       SetErrorCode("Block driver failed to make multi-transformation");
3176       return NULL;
3177     }
3178   }
3179   catch (Standard_Failure) {
3180     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3181     SetErrorCode(aFail->GetMessageString());
3182     return NULL;
3183   }
3184
3185   //Make a Python command
3186   GEOM::TPythonDump(aFunction) << aCopy << " = geompy.MakeMultiTransformation2D("
3187     << theObject << ", " << theDirFace1U << ", " << theDirFace2U << ", " << theNbTimesU
3188       << ", " << theDirFace1V << ", " << theDirFace2V << ", " << theNbTimesV << ")";
3189
3190   SetErrorCode(OK);
3191   return aCopy;
3192 }
3193
3194 //=============================================================================
3195 /*!
3196  *  Propagate
3197  */
3198 //=============================================================================
3199 Handle(TColStd_HSequenceOfTransient) GEOMImpl_IBlocksOperations::Propagate
3200                                                  (Handle(GEOM_Object) theShape)
3201 {
3202   SetErrorCode(KO);
3203
3204   if (theShape.IsNull()) return NULL;
3205
3206   TopoDS_Shape aShape = theShape->GetValue();
3207   if (aShape.IsNull()) return NULL;
3208
3209   TopTools_IndexedMapOfShape anIndices;
3210   TopExp::MapShapes(aShape, anIndices);
3211
3212   TopTools_IndexedDataMapOfShapeListOfShape MEW;
3213   GEOMImpl_Block6Explorer::MapShapesAndAncestors
3214     (aShape, TopAbs_EDGE, TopAbs_WIRE, MEW);
3215   Standard_Integer ie, nbEdges = MEW.Extent();
3216
3217   // Result
3218   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
3219
3220   TopTools_MapOfShape mapAcceptedEdges;
3221   TCollection_AsciiString aListRes, anEntry;
3222
3223   // Sort shapes in current chain (Mantis issue 21053)
3224   TopTools_DataMapOfShapeListOfShape aMapChains;
3225   TopTools_ListOfShape aFirstInChains;
3226
3227   for (ie = 1; ie <= nbEdges; ie++) {
3228     TopoDS_Shape curE = MEW.FindKey(ie);
3229
3230     if (mapAcceptedEdges.Contains(curE)) continue;
3231
3232     // Build the chain
3233     TopTools_ListOfShape currentChain;
3234     TopTools_ListOfShape listPrevEdges;
3235
3236     currentChain.Append(curE);
3237     listPrevEdges.Append(curE);
3238     mapAcceptedEdges.Add(curE);
3239
3240     // Collect all edges pass by pass
3241     while (listPrevEdges.Extent() > 0) {
3242       // List of edges, added to chain on this cycle pass
3243       TopTools_ListOfShape listCurEdges;
3244
3245       // Find the next portion of edges
3246       TopTools_ListIteratorOfListOfShape itE (listPrevEdges);
3247       for (; itE.More(); itE.Next()) {
3248         TopoDS_Shape anE = itE.Value();
3249
3250         // Iterate on faces, having edge <anE>
3251         TopTools_ListIteratorOfListOfShape itW (MEW.FindFromKey(anE));
3252         for (; itW.More(); itW.Next()) {
3253           TopoDS_Shape aW = itW.Value();
3254           TopoDS_Shape anOppE;
3255
3256           BRepTools_WireExplorer aWE (TopoDS::Wire(aW));
3257           Standard_Integer nb = 1, found = 0;
3258           TopTools_Array1OfShape anEdges (1,4);
3259           for (; aWE.More(); aWE.Next(), nb++) {
3260             if (nb > 4) {
3261               found = 0;
3262               break;
3263             }
3264             anEdges(nb) = aWE.Current();
3265             if (anEdges(nb).IsSame(anE)) found = nb;
3266           }
3267
3268           if (nb == 5 && found > 0) {
3269             // Quadrangle face found, get an opposite edge
3270             Standard_Integer opp = found + 2;
3271             if (opp > 4) opp -= 4;
3272             anOppE = anEdges(opp);
3273
3274             if (!mapAcceptedEdges.Contains(anOppE)) {
3275               // Add found edge to the chain
3276               currentChain.Append(anOppE);
3277               listCurEdges.Append(anOppE);
3278               mapAcceptedEdges.Add(anOppE);
3279             }
3280           } // if (nb == 5 && found > 0)
3281         } // for (; itF.More(); itF.Next())
3282       } // for (; itE.More(); itE.Next())
3283
3284       listPrevEdges = listCurEdges;
3285     } // while (listPrevEdges.Extent() > 0)
3286
3287     // Sort shapes in current chain (Mantis issue 21053)
3288     GEOMUtils::SortShapes(currentChain, Standard_False);
3289     aFirstInChains.Append(currentChain.First());
3290     aMapChains.Bind(currentChain.First(), currentChain);
3291   }
3292
3293   // Sort chains (Mantis issue 21053)
3294   GEOMUtils::SortShapes(aFirstInChains, Standard_False);
3295
3296   // Store sorted chains in the document
3297   TopTools_ListIteratorOfListOfShape aChainsIt (aFirstInChains);
3298   for (; aChainsIt.More(); aChainsIt.Next()) {
3299     TopoDS_Shape aFirstInChain = aChainsIt.Value();
3300     const TopTools_ListOfShape& currentChain = aMapChains.Find(aFirstInChain);
3301
3302     // Store the chain in the document
3303     Handle(TColStd_HArray1OfInteger) anArray =
3304       new TColStd_HArray1OfInteger (1, currentChain.Extent());
3305
3306     // Fill array of sub-shape indices
3307     TopTools_ListIteratorOfListOfShape itSub (currentChain);
3308     for (int index = 1; itSub.More(); itSub.Next(), ++index) {
3309       int id = anIndices.FindIndex(itSub.Value());
3310       anArray->SetValue(index, id);
3311     }
3312
3313     // Add a new group object
3314     Handle(GEOM_Object) aChain = GetEngine()->AddSubShape(theShape, anArray);
3315
3316     // Set a GROUP type
3317     aChain->SetType(GEOM_GROUP);
3318
3319     // Set a sub-shape type
3320     TDF_Label aFreeLabel = aChain->GetFreeLabel();
3321     TDataStd_Integer::Set(aFreeLabel, (Standard_Integer)TopAbs_EDGE);
3322
3323     // Add the chain to the result
3324     aSeq->Append(aChain);
3325
3326     //Make a Python command
3327     TDF_Tool::Entry(aChain->GetEntry(), anEntry);
3328     aListRes += anEntry + ", ";
3329   }
3330
3331   if (aSeq->IsEmpty()) {
3332     SetErrorCode("There are no quadrangle faces in the shape");
3333     return aSeq;
3334   }
3335
3336   aListRes.Trunc(aListRes.Length() - 2);
3337
3338   // The Propagation doesn't change object so no new function is required.
3339   Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
3340
3341   // Make a Python command
3342   GEOM::TPythonDump(aFunction, /*append=*/true)
3343     << "[" << aListRes.ToCString() << "] = geompy.Propagate(" << theShape << ")";
3344
3345   SetErrorCode(OK);
3346   return aSeq;
3347 }