src/AdvancedEngine/GEOMImpl_IAdvancedOperations.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 //  File   : GEOMImpl_IAdvancedOperations.cxx
  20 //  Author : Vadim SANDLER, Open CASCADE S.A.S. (vadim.sandler@opencascade.com)
  21
  22 #include "GEOMImpl_IAdvancedOperations.hxx"
  23
  24 #include <Basics_OCCTVersion.hxx>
  25
  26 #include <utilities.h>
  27 #include <OpUtil.hxx>
  28 #include <Utils_ExceptHandlers.hxx>
  29
  30 #include "GEOM_Function.hxx"
  31 #include "GEOM_PythonDump.hxx"
  32 #include "GEOMUtils.hxx"
  33 #include "GEOMAlgo_Splitter.hxx"
  34 #include "GEOMAlgo_FinderShapeOn1.hxx"
  35
  36 #include "GEOMImpl_Gen.hxx"
  37 #include "GEOMImpl_Types.hxx"
  38
  39 #include "GEOMImpl_IBasicOperations.hxx"
  40 #include "GEOMImpl_IBooleanOperations.hxx"
  41 #include "GEOMImpl_IShapesOperations.hxx"
  42 #include "GEOMImpl_ITransformOperations.hxx"
  43 #include "GEOMImpl_IBlocksOperations.hxx"
  44 #include "GEOMImpl_I3DPrimOperations.hxx"
  45 #include "GEOMImpl_ILocalOperations.hxx"
  46 #include "GEOMImpl_IHealingOperations.hxx"
  47 #include "GEOMImpl_IGroupOperations.hxx"
  48
  49 #include "GEOMImpl_GlueDriver.hxx"
  50 #include "GEOMImpl_PipeTShapeDriver.hxx"
  51 #include "GEOMImpl_IPipeTShape.hxx"
  52 #include "GEOMImpl_DividedDiskDriver.hxx"
  53 #include "GEOMImpl_IDividedDisk.hxx"
  54 #include <GEOMImpl_SmoothingSurfaceDriver.hxx>
  55 #include <GEOMImpl_ISmoothingSurface.hxx>
  56 /*@@ insert new functions before this line @@ do not remove this line @@ do not remove this line @@*/
  57
  58 #include <TDF_Tool.hxx>
  59 #include <TFunction_DriverTable.hxx>
  60 #include <TFunction_Driver.hxx>
  61 #include <TFunction_Logbook.hxx>
  62 #include <TNaming_CopyShape.hxx>
  63
  64 #include <TopExp.hxx>
  65 #include <TopExp_Explorer.hxx>
  66 #include <TopoDS.hxx>
  67 #include <TopoDS_Vertex.hxx>
  68 #include <TopTools_IndexedMapOfShape.hxx>
  69 #include <TopTools_ListIteratorOfListOfShape.hxx>
  70 #include <TColStd_IndexedDataMapOfTransientTransient.hxx>
  71
  72 #include <BRep_Builder.hxx>
  73 #include <BRep_Tool.hxx>
  74
  75 #include <BRepAdaptor_Surface.hxx>
  76 #include <BRepAlgoAPI_Cut.hxx>
  77 #include <BRepAlgoAPI_Fuse.hxx>
  78 #include <BRepBuilderAPI_MakeFace.hxx>
  79 #include <BRepBuilderAPI_MakeVertex.hxx>
  80 #include <BRepBuilderAPI_Transform.hxx>
  81 #include <BRepPrimAPI_MakeCone.hxx>
  82 #include <BRepPrimAPI_MakeCylinder.hxx>
  83
  84 #include <gp_Ax3.hxx>
  85 #include <gp_Pln.hxx>
  86 #include <gp_Pnt.hxx>
  87 #include <gp_Vec.hxx>
  88 #include <GC_MakeConicalSurface.hxx>
  89 #include <Geom_CylindricalSurface.hxx>
  90
  91 #include <ShapeAnalysis_Edge.hxx>
  92
  93 #include <cmath>
  94
  95 #include "AdvancedEngine_Types.hxx"
  96
  97 #include <Standard_Stream.hxx>
  98 #include <Standard_Failure.hxx>
  99 #include <StdFail_NotDone.hxx>
 100 #include <Standard_ErrorHandler.hxx> // CAREFUL ! position of this file is critic : see Lucien PIGNOLONI / OCC
 101
 102 #define HALF_LENGTH_MAIN_PIPE     "Main pipe half length" //"Tuyau principal - demi longueur"
 103 #define HALF_LENGTH_INCIDENT_PIPE "Incident pipe half length" //"Tuyau incident - demi longueur"
 104 #define CIRCULAR_QUARTER_PIPE     "Circular quarter of pipe" //"Circulaire - quart de tuyau"
 105 #define THICKNESS                 "Thickness" //"Epaisseur"
 106 #define FLANGE                    "Flange" // "Collerette"
 107 #define CHAMFER_OR_FILLET         "Chamfer or fillet" //"Chanfrein ou Raccord"
 108 #define JUNCTION_FACE_1           "Junction 1" //"Face de jonction 1"
 109 #define JUNCTION_FACE_2           "Junction 2" //"Face de jonction 2"
 110 #define JUNCTION_FACE_3           "Junction 3" //"Face de jonction 3"
 111
 112 #define FIND_GROUPS_BY_POINTS 1
 113
 114 //=============================================================================
 115 /*!
 116  *  Constructor
 117  */
 118 //=============================================================================
 119 GEOMImpl_IAdvancedOperations::GEOMImpl_IAdvancedOperations(GEOM_Engine* theEngine, int theDocID) :
 120   GEOM_IOperations(theEngine, theDocID)
 121 {
 122   MESSAGE("GEOMImpl_IAdvancedOperations::GEOMImpl_IAdvancedOperations");
 123   myBasicOperations     = new GEOMImpl_IBasicOperations(GetEngine(), GetDocID());
 124   myBooleanOperations   = new GEOMImpl_IBooleanOperations(GetEngine(), GetDocID());
 125   myShapesOperations    = new GEOMImpl_IShapesOperations(GetEngine(), GetDocID());
 126   myTransformOperations = new GEOMImpl_ITransformOperations(GetEngine(), GetDocID());
 127   myBlocksOperations    = new GEOMImpl_IBlocksOperations(GetEngine(), GetDocID());
 128   my3DPrimOperations    = new GEOMImpl_I3DPrimOperations(GetEngine(), GetDocID());
 129   myLocalOperations     = new GEOMImpl_ILocalOperations(GetEngine(), GetDocID());
 130   myHealingOperations   = new GEOMImpl_IHealingOperations(GetEngine(), GetDocID());
 131   myGroupOperations     = new GEOMImpl_IGroupOperations(GetEngine(), GetDocID());
 132 }
 133
 134 //=============================================================================
 135 /*!
 136  *  Destructor
 137  */
 138 //=============================================================================
 139 GEOMImpl_IAdvancedOperations::~GEOMImpl_IAdvancedOperations()
 140 {
 141   MESSAGE("GEOMImpl_IAdvancedOperations::~GEOMImpl_IAdvancedOperations");
 142   delete myBasicOperations;
 143   delete myBooleanOperations;
 144   delete myShapesOperations;
 145   delete myTransformOperations;
 146   delete myBlocksOperations;
 147   delete my3DPrimOperations;
 148   delete myLocalOperations;
 149   delete myHealingOperations;
 150   delete myGroupOperations;
 151 }
 152
 153 //=============================================================================
 154 /*!
 155  *  SetPosition
 156  */
 157 //=============================================================================
 158 gp_Trsf GEOMImpl_IAdvancedOperations::GetPositionTrsf(double theL1, double theL2,
 159                                                       Handle(GEOM_Object) theP1,
 160                                                       Handle(GEOM_Object) theP2,
 161                                                       Handle(GEOM_Object) theP3)
 162 {
 163   // Old Local Coordinates System oldLCS
 164   gp_Pnt P0(0, 0, 0);
 165   gp_Pnt P1(-theL1, 0, 0);
 166   gp_Pnt P2(theL1, 0, 0);
 167   gp_Pnt P3(0, 0, theL2);
 168
 169   gp_Dir oldX(gp_Vec(P1, P2));
 170   gp_Dir oldZ(gp_Vec(P0, P3));
 171   gp_Ax3 oldLCS(P0, oldZ, oldX);
 172
 173   // New Local Coordinates System newLCS
 174   double LocX, LocY, LocZ;
 175   gp_Pnt newP1 = BRep_Tool::Pnt(TopoDS::Vertex(theP1->GetValue()));
 176   gp_Pnt newP2 = BRep_Tool::Pnt(TopoDS::Vertex(theP2->GetValue()));
 177   gp_Pnt newP3 = BRep_Tool::Pnt(TopoDS::Vertex(theP3->GetValue()));
 178   LocX = (newP1.X() + newP2.X()) / 2.;
 179   LocY = (newP1.Y() + newP2.Y()) / 2.;
 180   LocZ = (newP1.Z() + newP2.Z()) / 2.;
 181   gp_Pnt newO(LocX, LocY, LocZ);
 182
 183   gp_Dir newX(gp_Vec(newP1, newP2)); // P1P2 Vector
 184   gp_Dir newZ(gp_Vec(newO, newP3)); // OP3 Vector
 185   gp_Ax3 newLCS = gp_Ax3(newO, newZ, newX);
 186
 187   gp_Trsf aTrsf;
 188   aTrsf.SetDisplacement(oldLCS, newLCS);
 189
 190   return aTrsf;
 191 }
 192
 193 //=============================================================================
 194 /*!
 195  *  CheckCompatiblePosition
 196  *
 197  */
 198 //=============================================================================
 199 bool GEOMImpl_IAdvancedOperations::CheckCompatiblePosition(double& theL1, double& theL2,
 200                                                            Handle(GEOM_Object) theP1,
 201                                                            Handle(GEOM_Object) theP2,
 202                                                            Handle(GEOM_Object) theP3,
 203                                                            double theTolerance)
 204 {
 205   SetErrorCode(KO);
 206   gp_Pnt P1 = BRep_Tool::Pnt(TopoDS::Vertex(theP1->GetValue()));
 207   gp_Pnt P2 = BRep_Tool::Pnt(TopoDS::Vertex(theP2->GetValue()));
 208   gp_Pnt P3 = BRep_Tool::Pnt(TopoDS::Vertex(theP3->GetValue()));
 209
 210   double d12 = P1.Distance(P2);
 211   double d13 = P1.Distance(P3);
 212   double d23 = P2.Distance(P3);
 213   //    double d2 = newO.Distance(P3);
 214
 215   if (Abs(d12) <= Precision::Confusion()) {
 216     SetErrorCode("Junctions points P1 and P2 are identical");
 217     return false;
 218   }
 219   if (Abs(d13) <= Precision::Confusion()) {
 220     SetErrorCode("Junctions points P1 and P3 are identical");
 221     return false;
 222   }
 223   if (Abs(d23) <= Precision::Confusion()) {
 224     SetErrorCode("Junctions points P2 and P3 are identical");
 225     return false;
 226   }
 227
 228
 229   double newL1 = 0.5 * d12;
 230   double newL2 = sqrt(pow(d13,2)-pow(newL1,2));
 231   //
 232   // theL1*(1-theTolerance) <= newL1 <= theL1*(1+theTolerance)
 233   //
 234   if (fabs(newL1 - theL1) > Precision::Approximation()) {
 235     if ( (newL1 * (1 - theTolerance) -theL1 <= Precision::Approximation()) &&
 236          (newL1 * (1 + theTolerance) -theL1 >= Precision::Approximation()) ) {
 237       //            std::cerr << "theL1 = newL1" << std::endl;
 238       theL1 = newL1;
 239     } else {
 240       theL1 = -1;
 241       SetErrorCode("Dimension for main pipe (L1) is incompatible with new position");
 242       return false;
 243     }
 244   }
 245
 246   //
 247   // theL2*(1-theTolerance) <= newL2  <= theL2*(1+theTolerance)
 248   //
 249   if (fabs(newL2 - theL2) > Precision::Approximation()) {
 250     if ( (newL2 * (1 - theTolerance) -theL2 <= Precision::Approximation()) &&
 251          (newL2 * (1 + theTolerance) -theL2 >= Precision::Approximation()) ) {
 252       theL2 = newL2;
 253     } else {
 254       theL2 = -1;
 255       SetErrorCode("Dimension for incident pipe (L2) is incompatible with new position");
 256       return false;
 257     }
 258   }
 259
 260   SetErrorCode(OK);
 261   return true;
 262
 263 }
 264
 265 //=============================================================================
 266 /*!
 267  *  Generate the propagation groups of a Pipe T-Shape used for hexa mesh
 268  */
 269 //=============================================================================
 270 bool GEOMImpl_IAdvancedOperations::MakeGroups(Handle(GEOM_Object) theShape, int shapeType,
 271                                               double theR1, double theW1, double theL1,
 272                                               double theR2, double theW2, double theL2,
 273                                               double theH, double theW, double theRF,
 274                                               Handle(TColStd_HSequenceOfTransient) theSeq,
 275                                               gp_Trsf aTrsf)
 276 {
 277   SetErrorCode(KO);
 278
 279   if (theShape.IsNull()) return false;
 280
 281   TopoDS_Shape aShape = theShape->GetValue();
 282   if (aShape.IsNull()) {
 283     SetErrorCode("Shape is not defined");
 284     return false;
 285   }
 286
 287   gp_Trsf aTrsfInv = aTrsf.Inverted();
 288
 289 //   int expectedGroups = 0;
 290 //   if (shapeType == TSHAPE_BASIC)
 291 //     if (Abs(theR2+theW2-theR1-theW1) <= Precision::Approximation())
 292 //       expectedGroups = 10;
 293 //     else
 294 //       expectedGroups = 11;
 295 //   else if (shapeType == TSHAPE_CHAMFER || shapeType == TSHAPE_FILLET)
 296 //     expectedGroups = 12;
 297
 298   double aR1Ext = theR1 + theW1;
 299   double aR2Ext = theR2 + theW2;
 300
 301   /////////////////////////
 302   //// Groups of Faces ////
 303   /////////////////////////
 304
 305   //
 306   // Comment the following lines when GetInPlace bug is solved
 307   // == BEGIN
 308   // Workaround of GetInPlace bug
 309   // Create a bounding box that fits the shape
 310   Handle(GEOM_Object) aBox = my3DPrimOperations->MakeBoxDXDYDZ(2*theL1, 2*aR1Ext, aR1Ext+theL2);
 311   aBox->GetLastFunction()->SetDescription("");
 312   myTransformOperations->TranslateDXDYDZ(aBox, -theL1, -aR1Ext, -aR1Ext);
 313   aBox->GetLastFunction()->SetDescription("");
 314   // Apply transformation to box
 315   BRepBuilderAPI_Transform aTransformationBox(aBox->GetValue(), aTrsf, Standard_False);
 316   TopoDS_Shape aBoxShapeTrsf = aTransformationBox.Shape();
 317   aBox->GetLastFunction()->SetValue(aBoxShapeTrsf);
 318
 319   // Get the shell of the box
 320   Handle(GEOM_Object) aShell = Handle(GEOM_Object)::DownCast
 321     (myShapesOperations->MakeExplode(aBox, TopAbs_SHELL, true)->Value(1));
 322   aBox->GetLastFunction()->SetDescription("");
 323   aShell->GetLastFunction()->SetDescription("");
 324   // Get the common shapes between shell and shape
 325   Handle(GEOM_Object) aCommonCompound = myBooleanOperations->MakeBoolean
 326                             (theShape, aShell, 1, Standard_False); // MakeCommon
 327   if (aCommonCompound.IsNull()) {
 328     SetErrorCode(myBooleanOperations->GetErrorCode());
 329     return false;
 330   }
 331   aCommonCompound->GetLastFunction()->SetDescription("");
 332   // Explode the faces of common shapes => 3 faces
 333   Handle(TColStd_HSequenceOfTransient) aCommonFaces =
 334     myShapesOperations->MakeExplode(aCommonCompound, TopAbs_FACE, true);
 335   aCommonCompound->GetLastFunction()->SetDescription("");
 336   std::list<Handle(GEOM_Object)> aCompoundOfFacesList;
 337
 338   for (int i=0 ; i<= aCommonFaces->Length()-4 ; i+=4) {
 339     std::list<Handle(GEOM_Object)> aFacesList;
 340     for (int j = 1 ; j <= 4 ; j++) {
 341       Handle(GEOM_Object) aFace = Handle(GEOM_Object)::DownCast(aCommonFaces->Value(i+j)); // Junction faces
 342       if (!aFace.IsNull()) {
 343         aFace->GetLastFunction()->SetDescription("");
 344         aFacesList.push_back(aFace);
 345       }
 346     }
 347     Handle(GEOM_Object) aCompoundOfFaces = myShapesOperations->MakeCompound(aFacesList);
 348     if (!aCompoundOfFaces.IsNull()) {
 349       aCompoundOfFaces->GetLastFunction()->SetDescription("");
 350       aCompoundOfFacesList.push_back(aCompoundOfFaces);
 351     }
 352   }
 353
 354   if (aCompoundOfFacesList.size() == 3) {
 355     Handle(GEOM_Object) aPln1 = aCompoundOfFacesList.front();
 356     aCompoundOfFacesList.pop_front();
 357     Handle(GEOM_Object) aPln2 = aCompoundOfFacesList.front();
 358     aCompoundOfFacesList.pop_front();
 359     Handle(GEOM_Object) aPln3 = aCompoundOfFacesList.front();
 360     aCompoundOfFacesList.pop_front();
 361     // == END
 362     //
 363
 364
 365     //     Uncomment the following lines when GetInPlace bug is solved
 366     //     == BEGIN
 367 //     Handle(GEOM_Object) aP1 = myBasicOperations->MakePointXYZ(-theL1, 0, 0);
 368 //     Handle(GEOM_Object) aP2 = myBasicOperations->MakePointXYZ(-0, 0, theL2);
 369 //     Handle(GEOM_Object) aP3 = myBasicOperations->MakePointXYZ(theL1, 0, 0);
 370 //     aP1->GetLastFunction()->SetDescription("");
 371 //     aP2->GetLastFunction()->SetDescription("");
 372 //     aP3->GetLastFunction()->SetDescription("");
 373 //     Handle(GEOM_Object) aV1 = myBasicOperations->MakeVectorDXDYDZ(-1, 0, 0);
 374 //     Handle(GEOM_Object) aV2 = myBasicOperations->MakeVectorDXDYDZ(0, 0, 1);
 375 //     Handle(GEOM_Object) aV3 = myBasicOperations->MakeVectorDXDYDZ(1, 0, 0);
 376 //     aV1->GetLastFunction()->SetDescription("");
 377 //     aV2->GetLastFunction()->SetDescription("");
 378 //     aV3->GetLastFunction()->SetDescription("");
 379 //     Handle(GEOM_Object) aPln1 = myBasicOperations->MakePlanePntVec(aP1, aV1, 2*(aR1Ext+theL2));
 380 //     Handle(GEOM_Object) aPln2 = myBasicOperations->MakePlanePntVec(aP2, aV2, 2*(aR2Ext));
 381 //     Handle(GEOM_Object) aPln3 = myBasicOperations->MakePlanePntVec(aP3, aV3, 2*(aR1Ext+theL2));
 382 //     aPln1->GetLastFunction()->SetDescription("");
 383 //     aPln2->GetLastFunction()->SetDescription("");
 384 //     aPln3->GetLastFunction()->SetDescription("");
 385 //
 386 //     BRepBuilderAPI_Transform aTransformation1(aPln1->GetValue(), aTrsf, Standard_False);
 387 //     TopoDS_Shape aTrsf_Shape1 = aTransformation1.Shape();
 388 //     aPln1->GetLastFunction()->SetValue(aTrsf_Shape1);
 389 //     BRepBuilderAPI_Transform aTransformation2(aPln2->GetValue(), aTrsf, Standard_False);
 390 //     TopoDS_Shape aTrsf_Shape2 = aTransformation2.Shape();
 391 //     aPln2->GetLastFunction()->SetValue(aTrsf_Shape2);
 392 //     BRepBuilderAPI_Transform aTransformation3(aPln3->GetValue(), aTrsf, Standard_False);
 393 //     TopoDS_Shape aTrsf_Shape3 = aTransformation3.Shape();
 394 //     aPln3->GetLastFunction()->SetValue(aTrsf_Shape3);
 395     //     == END
 396     //
 397
 398     Handle(GEOM_Object) junctionFaces1 = myShapesOperations->GetInPlace(theShape, aPln1);
 399     if (junctionFaces1.IsNull())
 400       junctionFaces1 = myShapesOperations->GetShapesOnShapeAsCompound
 401         (aPln1, theShape, TopAbs_FACE,  GEOMAlgo_ST_ONIN);
 402     if (!junctionFaces1.IsNull()) {
 403       junctionFaces1->GetLastFunction()->SetDescription("");
 404       junctionFaces1->SetName("JUNCTION_FACE_1");
 405       theSeq->Append(junctionFaces1);
 406     }
 407     else {
 408       SetErrorCode("Junction face 1 not found");
 409       //        theSeq->Append(aPln1);
 410       //        return false;
 411     }
 412     Handle(GEOM_Object) junctionFaces2 = myShapesOperations->GetInPlace(theShape, aPln2);
 413     if (junctionFaces2.IsNull())
 414       junctionFaces2 = myShapesOperations->GetShapesOnShapeAsCompound
 415         (aPln2, theShape, TopAbs_FACE,  GEOMAlgo_ST_ONIN);
 416     if (!junctionFaces2.IsNull()) {
 417       junctionFaces2->GetLastFunction()->SetDescription("");
 418       junctionFaces2->SetName("JUNCTION_FACE_2");
 419       theSeq->Append(junctionFaces2);
 420     }
 421     else {
 422       SetErrorCode("Junction face 2 not found");
 423       //        theSeq->Append(aPln2);
 424       //        return false;
 425     }
 426     Handle(GEOM_Object) junctionFaces3 = myShapesOperations->GetInPlace(theShape, aPln3);
 427     if (junctionFaces3.IsNull())
 428       junctionFaces3 = myShapesOperations->GetShapesOnShapeAsCompound
 429         (aPln3, theShape, TopAbs_FACE,  GEOMAlgo_ST_ONIN);
 430     if (!junctionFaces3.IsNull()) {
 431       junctionFaces3->GetLastFunction()->SetDescription("");
 432       junctionFaces3->SetName("JUNCTION_FACE_3");
 433       theSeq->Append(junctionFaces3);
 434     }
 435     else {
 436       SetErrorCode("Junction face 3 not found");
 437       //        theSeq->Append(aPln3);
 438       //        return false;
 439     }
 440   // Comment the following lines when GetInPlace bug is solved
 441   // == BEGIN
 442   }
 443   //     == END
 444
 445   /////////////////////////
 446   //// Groups of Edges ////
 447   /////////////////////////
 448   // Result of propagate
 449
 450   Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
 451
 452   TCollection_AsciiString theDesc = aFunction->GetDescription();
 453   Handle(TColStd_HSequenceOfTransient) aSeqPropagate = myBlocksOperations->Propagate(theShape);
 454   if (aSeqPropagate.IsNull() || aSeqPropagate->Length() == 0) {
 455     SetErrorCode("Propagation groups not found");
 456     return false;
 457   }
 458   Standard_Integer aNbGroups = aSeqPropagate->Length();
 459   // Recover previous description to get rid of Propagate dump
 460   aFunction->SetDescription(theDesc);
 461
 462 #ifdef FIND_GROUPS_BY_POINTS
 463   // BEGIN: new groups search
 464
 465   //              W2  R2
 466   //            .----.-----.----.
 467   //           e|    |  |  |    |
 468   //            |    |  |  |    |
 469   //            .    |  |  |    .
 470   //         g / ''..|  |  |..'' \
 471   //       f  /      '''''''      \
 472   //  .---.--'..     |  |  |     ..'--.---.
 473   //  |a    \   '''...........'''   /     |
 474   //  |-------\------'  |  '------/-------.
 475   //  |         \       |       /         |
 476   // c|           \     |     /           |
 477   //  |    R1       \   |   /             |
 478   //  |               \ | /               |
 479   //  ._________________|_________________.
 480   //  |       L1        |                 |
 481   //  |                 |                 |
 482   //  |                 |                 |
 483   // b|                 |                 |
 484   //  |                 |                 |
 485   //  |-----------------|-----------------|
 486   //  |    W1           |                 |
 487   //  '-----------------'-----------------'
 488   //          d
 489
 490   // "Thickness" group (a)
 491   gp_Pnt aPntA (-theL1, 0, theR1 + theW1/2.);
 492   aPntA.Transform(aTrsf);
 493   BRepBuilderAPI_MakeVertex mkVertexA (aPntA);
 494   TopoDS_Vertex aVertA = TopoDS::Vertex(mkVertexA.Shape());
 495   TopoDS_Shape anEdgeA = GEOMUtils::GetEdgeNearPoint(aShape, aVertA);
 496
 497   // "Circular quarter of pipe" group (b)
 498   gp_Pnt aPntB (-theL1, -aR1Ext * sin(M_PI/4.), -aR1Ext * sin(M_PI/4.));
 499   aPntB.Transform(aTrsf);
 500   BRepBuilderAPI_MakeVertex mkVertexB (aPntB);
 501   TopoDS_Vertex aVertB = TopoDS::Vertex(mkVertexB.Shape());
 502   TopoDS_Shape anEdgeB = GEOMUtils::GetEdgeNearPoint(aShape, aVertB);
 503
 504   // "Circular quarter of pipe" group (c)
 505   gp_Pnt aPntC (-theL1, -aR1Ext * sin(M_PI/4.), aR1Ext * sin(M_PI/4.));
 506   aPntC.Transform(aTrsf);
 507   BRepBuilderAPI_MakeVertex mkVertexC (aPntC);
 508   TopoDS_Vertex aVertC = TopoDS::Vertex(mkVertexC.Shape());
 509   TopoDS_Shape anEdgeC = GEOMUtils::GetEdgeNearPoint(aShape, aVertC);
 510
 511   // "Main pipe half length" group (d)
 512   gp_Pnt aPntD (-theL1/2., 0, -aR1Ext);
 513   aPntD.Transform(aTrsf);
 514   BRepBuilderAPI_MakeVertex mkVertexD (aPntD);
 515   TopoDS_Vertex aVertD = TopoDS::Vertex(mkVertexD.Shape());
 516   TopoDS_Shape anEdgeD = GEOMUtils::GetEdgeNearPoint(aShape, aVertD);
 517
 518   // "Incident pipe half length" group (e)
 519   double aTol10 = Precision::Confusion() * 10.;
 520   gp_Pnt aPntE (-aR2Ext, 0, theL2 - aTol10);
 521   aPntE.Transform(aTrsf);
 522   BRepBuilderAPI_MakeVertex mkVertexE (aPntE);
 523   TopoDS_Vertex aVertE = TopoDS::Vertex(mkVertexE.Shape());
 524   TopoDS_Shape anEdgeE = GEOMUtils::GetEdgeNearPoint(aShape, aVertE);
 525
 526   // "Flange" group (f)
 527   double aFx = - aR2Ext - aTol10;
 528   if (shapeType == TSHAPE_CHAMFER)
 529     aFx -= theW;
 530   else if (shapeType == TSHAPE_FILLET)
 531     aFx -= theRF;
 532   gp_Pnt aPntF (aFx, 0, aR1Ext);
 533   aPntF.Transform(aTrsf);
 534   BRepBuilderAPI_MakeVertex mkVertexF (aPntF);
 535   TopoDS_Vertex aVertF = TopoDS::Vertex(mkVertexF.Shape());
 536   TopoDS_Shape anEdgeF = GEOMUtils::GetEdgeNearPoint(aShape, aVertF);
 537
 538   // "Chamfer or Fillet" group (g)
 539   TopoDS_Shape anEdgeG;
 540   if (shapeType == TSHAPE_CHAMFER) {
 541     gp_Pnt aPntG (-aR2Ext - theW/2., 0, aR1Ext + theH/2.);
 542     aPntG.Transform(aTrsf);
 543     BRepBuilderAPI_MakeVertex mkVertexG (aPntG);
 544     TopoDS_Vertex aVertG = TopoDS::Vertex(mkVertexG.Shape());
 545     anEdgeG = GEOMUtils::GetEdgeNearPoint(aShape, aVertG);
 546   }
 547   else if (shapeType == TSHAPE_FILLET) {
 548     gp_Pnt aPntG (-aR2Ext - theRF/2., 0, aR1Ext + theRF/2.);
 549     aPntG.Transform(aTrsf);
 550     BRepBuilderAPI_MakeVertex mkVertexG (aPntG);
 551     TopoDS_Vertex aVertG = TopoDS::Vertex(mkVertexG.Shape());
 552     anEdgeG = GEOMUtils::GetEdgeNearPoint(aShape, aVertG);
 553   }
 554
 555   for (int i = 1 ; i <= aNbGroups; i++) {
 556     Handle(GEOM_Object) aGroup = Handle(GEOM_Object)::DownCast(aSeqPropagate->Value(i));
 557     if (aGroup.IsNull())
 558       continue;
 559
 560     TopoDS_Shape aGroupShape = aGroup->GetValue();
 561     TopTools_IndexedMapOfShape anEdgesMap;
 562     TopExp::MapShapes(aGroupShape, TopAbs_EDGE, anEdgesMap);
 563
 564     if (anEdgesMap.Contains(anEdgeA)) { // a
 565       aGroup->SetName("THICKNESS");
 566       theSeq->Append(aGroup);
 567     }
 568     else if (anEdgesMap.Contains(anEdgeB)) { // b
 569       aGroup->SetName("CIRCULAR_QUARTER_PIPE");
 570       theSeq->Append(aGroup);
 571     }
 572     else if (anEdgesMap.Contains(anEdgeC)) { // c
 573       aGroup->SetName("CIRCULAR_QUARTER_PIPE");
 574       theSeq->Append(aGroup);
 575     }
 576     else if (anEdgesMap.Contains(anEdgeD)) { // d
 577       aGroup->SetName("HALF_LENGTH_MAIN_PIPE");
 578       theSeq->Append(aGroup);
 579     }
 580     else if (anEdgesMap.Contains(anEdgeE)) { // e
 581       aGroup->SetName("HALF_LENGTH_INCIDENT_PIPE");
 582       theSeq->Append(aGroup);
 583     }
 584     else if (anEdgesMap.Contains(anEdgeF)) { // f
 585       aGroup->SetName("FLANGE");
 586       theSeq->Append(aGroup);
 587     }
 588     else if (shapeType == TSHAPE_CHAMFER) { // g
 589       if (anEdgesMap.Contains(anEdgeG)) {
 590         aGroup->SetName("CHAMFER");
 591         theSeq->Append(aGroup);
 592       }
 593     }
 594     else if (shapeType == TSHAPE_FILLET) { // g
 595       if (anEdgesMap.Contains(anEdgeG)) {
 596         aGroup->SetName("FILLET");
 597         theSeq->Append(aGroup);
 598       }
 599     }
 600     else {
 601     }
 602   }
 603   // END: new groups search
 604 #else
 605   bool addGroup;
 606   bool circularFoundAndAdded = false;
 607   bool circularFound10 = false;
 608   bool incidentPipeFound = false;
 609   bool mainPipeFound = false;
 610   bool mainPipeFoundAndAdded = false;
 611   bool radialFound =false;
 612   bool flangeFound = false;
 613   bool flangeFoundAndAdded = false;
 614   bool chamferOrFilletFound = false;
 615
 616   for (int i = 1 ; i <= aNbGroups; i++) {
 617     addGroup = false;
 618
 619     Handle(GEOM_Object) aGroup = Handle(GEOM_Object)::DownCast(aSeqPropagate->Value(i));
 620     if (aGroup.IsNull())
 621       continue;
 622
 623     TopoDS_Shape aGroupShape = aGroup->GetValue();
 624     BRepBuilderAPI_Transform aTransformationShapeInv (aGroupShape, aTrsfInv, Standard_False);
 625     TopoDS_Shape aGroupShapeTrsfInv = aTransformationShapeInv.Shape();
 626
 627     TopTools_IndexedMapOfShape anEdgesMap;
 628     TopExp::MapShapes(aGroupShapeTrsfInv,TopAbs_EDGE, anEdgesMap);
 629     Standard_Integer nbEdges = anEdgesMap.Extent();
 630
 631     if (shapeType == TSHAPE_BASIC) {
 632       if ((nbEdges >= 21) || /*R1Ext = R2Ext*/(nbEdges == 17)) { // 17, 17+8*{1,2,3}, 21, 21+8*{1,2,3}
 633         addGroup = true;
 634         aGroup->SetName("THICKNESS");
 635       }
 636       else if (nbEdges == 6) {
 637         if (!circularFoundAndAdded) {
 638           circularFoundAndAdded = true;
 639           addGroup = true;
 640           aGroup->SetName("CIRCULAR_QUARTER_PIPE");
 641         }
 642       }
 643       else if (nbEdges == 8) {
 644         incidentPipeFound = true;
 645         mainPipeFound = false;
 646         radialFound = false;
 647         flangeFound = false;
 648
 649         TopExp_Explorer Ex(aGroupShapeTrsfInv,TopAbs_VERTEX);
 650         while (Ex.More()) {
 651           gp_Pnt aP =  BRep_Tool::Pnt(TopoDS::Vertex(Ex.Current()));
 652           double x=aP.X(), y=aP.Y(), z=aP.Z();
 653
 654
 655           if ((Abs(x) > aR2Ext + Precision::Confusion()) ||
 656               (Abs(y) > aR2Ext + Precision::Confusion())) {
 657             incidentPipeFound = false;
 658           }
 659
 660           if ( z < -Precision::Confusion()) {
 661             // length of main pipe
 662             mainPipeFound = true;
 663             if (!mainPipeFoundAndAdded) {
 664               mainPipeFoundAndAdded = true;
 665               addGroup = true;
 666               aGroup->SetName("HALF_LENGTH_MAIN_PIPE");
 667             }
 668           }
 669
 670           else if (Abs(x) > (theL1-Precision::Confusion())) {
 671             // discretisation circulaire
 672             radialFound = true;
 673             if (!circularFoundAndAdded) {
 674               circularFoundAndAdded = true;
 675               addGroup = true;
 676               aGroup->SetName("CIRCULAR_QUARTER_PIPE");
 677             }
 678           }
 679           Ex.Next();
 680         }
 681         if (incidentPipeFound) {
 682           addGroup = true;
 683           aGroup->SetName("HALF_LENGTH_INCIDENT_PIPE");
 684         }
 685         if (!addGroup && (!incidentPipeFound &&
 686                           !radialFound &&
 687                           !mainPipeFound &&
 688                           !flangeFound)) {
 689           // Flange (collerette)
 690           flangeFound = true;
 691           addGroup = true;
 692           aGroup->SetName("FLANGE");
 693         }
 694       }
 695       else
 696         continue;
 697     }
 698     else if (shapeType == TSHAPE_CHAMFER || shapeType == TSHAPE_FILLET) {
 699       if (nbEdges >= 25) { // 25, 25+8, 25+16, 25+24
 700         addGroup = true;
 701         aGroup->SetName("THICKNESS");
 702       }
 703       else if ((nbEdges == 10) || (nbEdges == 6)) {
 704         if (!circularFoundAndAdded) {
 705           addGroup = true;
 706           circularFoundAndAdded = true;
 707           aGroup->SetName("CIRCULAR_QUARTER_PIPE");
 708           if (nbEdges == 10) {
 709             circularFound10 = true;
 710           }
 711         }
 712         else if (!circularFound10 && nbEdges == 10) {
 713           circularFound10 = true;
 714           addGroup = true;
 715           aGroup->SetName("CIRCULAR_QUARTER_PIPE");
 716         }
 717       }
 718       else if (nbEdges == 8) {
 719         incidentPipeFound = true;
 720         mainPipeFound = true;
 721         flangeFound = false;
 722
 723         bool isNearZ0 = false;
 724         bool isBelowZ0 = false;
 725
 726         TopExp_Explorer Ex (aGroupShapeTrsfInv,TopAbs_VERTEX);
 727         while (Ex.More()) {
 728           gp_Pnt aP = BRep_Tool::Pnt(TopoDS::Vertex(Ex.Current()));
 729           double x=aP.X(), y=aP.Y(), z=aP.Z();
 730
 731           // tuy_princ_long_avant & tuy_princ_long_apres
 732           //bool isMain = (((z < Precision::Confusion()) || (x < Precision::Confusion())) &&
 733           //               ((y <= aR1Ext + Precision::Confusion()) ||
 734           //                (y <= -(aR1Ext + Precision::Confusion())) ||
 735           //                (y <= theR1 + Precision::Confusion()) ||
 736           //                (y == -(theR1 + Precision::Confusion()))));
 737           bool isMain = ((z < Precision::Confusion() || x < Precision::Confusion()) &&
 738                          (fabs(y) > theR1 - Precision::Confusion() ||
 739                           fabs(y) < Precision::Confusion()));
 740
 741           if (!isMain) {
 742             mainPipeFound = false;
 743           }
 744
 745           // collerette
 746           //if (z < Precision::Confusion() && !isMain) {
 747           //  flangeFound = true;
 748           //  if (!flangeFoundAndAdded) {
 749           //    flangeFoundAndAdded = true;
 750           //    addGroup = true;
 751           //    aGroup->SetName("FLANGE");
 752           //  }
 753           //}
 754           if (fabs(z) < Precision::Confusion()) isNearZ0 = true;
 755           if (z < - Precision::Confusion()) isBelowZ0 = true;
 756
 757           // tuyau incident
 758           if ((Abs(x) > aR2Ext + Precision::Confusion()) ||
 759               (Abs(y) > aR2Ext + Precision::Confusion())) {
 760             incidentPipeFound = false;
 761           }
 762           Ex.Next();
 763         }
 764         if (mainPipeFound) {
 765           addGroup = true;
 766           aGroup->SetName("HALF_LENGTH_MAIN_PIPE");
 767         }
 768         if (incidentPipeFound) {
 769           addGroup = true;
 770           aGroup->SetName("HALF_LENGTH_INCIDENT_PIPE");
 771         }
 772         if (isNearZ0 && !isBelowZ0) {
 773           flangeFound = true;
 774           if (!flangeFoundAndAdded) {
 775             flangeFoundAndAdded = true;
 776             addGroup = true;
 777             aGroup->SetName("FLANGE");
 778           }
 779         }
 780         if (!addGroup && (!incidentPipeFound &&
 781                           !mainPipeFound &&
 782                           !flangeFound &&
 783                           !chamferOrFilletFound)) {
 784           addGroup = true;
 785           chamferOrFilletFound = true;
 786           if (shapeType == TSHAPE_CHAMFER)
 787             aGroup->SetName("CHAMFER");
 788           else
 789             aGroup->SetName("FILLET");
 790         }
 791       }
 792       else
 793         continue;
 794     }
 795     // Add group to the list
 796     if (addGroup)
 797       theSeq->Append(aGroup);
 798   }
 799 #endif
 800
 801   SetErrorCode(OK);
 802   return true;
 803 }
 804
 805 //=============================================================================
 806 /*!
 807  *  Return faces that are laying on surface.
 808  */
 809 //=============================================================================
 810 bool GEOMImpl_IAdvancedOperations::GetFacesOnSurf
 811                      (const TopoDS_Shape &theShape,
 812                       const Handle_Geom_Surface& theSurface,
 813                       const Standard_Real theTolerance,
 814                       TopTools_ListOfShape &theFaces)
 815 {
 816   GEOMAlgo_FinderShapeOn1 aFinder;
 817
 818   aFinder.SetShape(theShape);
 819   aFinder.SetTolerance(theTolerance);
 820   aFinder.SetSurface(theSurface);
 821   aFinder.SetShapeType(TopAbs_FACE);
 822   aFinder.SetState(GEOMAlgo_ST_ON);
 823
 824   // Sets the minimal number of inner points for the faces that do not have own
 825   // inner points at all (for e.g. rectangular planar faces have just 2 triangles).
 826   // Default value=3
 827   aFinder.SetNbPntsMin(3);
 828   // Sets the maximal number of inner points for edges or faces.
 829   // It is usefull for the cases when this number is very big (e.g =2000) to improve
 830   // the performance. If this value =0, all inner points will be taken into account.
 831   // Default value=0
 832   aFinder.SetNbPntsMax(100);
 833   aFinder.Perform();
 834
 835   // Interprete results
 836   Standard_Integer iErr = aFinder.ErrorStatus();
 837   // the detailed description of error codes is in GEOMAlgo_FinderShapeOn1.cxx
 838   if (iErr) {
 839     MESSAGE(" iErr : " << iErr);
 840     TCollection_AsciiString aMsg (" iErr : ");
 841     aMsg += TCollection_AsciiString(iErr);
 842     SetErrorCode(aMsg);
 843     return false;
 844   }
 845   Standard_Integer iWrn = aFinder.WarningStatus();
 846   // the detailed description of warning codes is in GEOMAlgo_FinderShapeOn1.cxx
 847   if (iWrn) {
 848     MESSAGE(" *** iWrn : " << iWrn);
 849   }
 850
 851   const TopTools_ListOfShape &aListRes = aFinder.Shapes(); // the result
 852   TopTools_ListIteratorOfListOfShape anIter (aListRes);
 853
 854   for (; anIter.More(); anIter.Next()) {
 855     theFaces.Append(anIter.Value());
 856   }
 857
 858   return true;
 859 }
 860
 861 //=============================================================================
 862 /*!
 863  *  Creates and returns conical face.
 864  */
 865 //=============================================================================
 866 TopoDS_Shape GEOMImpl_IAdvancedOperations::MakeConicalFace
 867                                   (const gp_Ax2 &theAxis,
 868                                    const double theRadius,
 869                                    const double theRadiusThin,
 870                                    const double theHeight,
 871                                    const gp_Trsf &theTrsf)
 872 {
 873   BRepPrimAPI_MakeCone aMkCone (theAxis, theRadius, theRadiusThin, theHeight);
 874   TopoDS_Shape aResult;
 875
 876   aMkCone.Build();
 877   if (aMkCone.IsDone()) {
 878     TopExp_Explorer anExp(aMkCone.Shape(), TopAbs_FACE);
 879
 880     for (; anExp.More(); anExp.Next()) {
 881       TopoDS_Face aFace = TopoDS::Face(anExp.Current());
 882
 883       if (aFace.IsNull() == Standard_False) {
 884         BRepAdaptor_Surface anAdaptor(aFace, Standard_False);
 885
 886         if (anAdaptor.GetType() == GeomAbs_Cone) {
 887           // This is a conical face. Transform and return it.
 888           BRepBuilderAPI_Transform aTransf(aFace, theTrsf, Standard_False);
 889
 890           aResult = aTransf.Shape();
 891           break;
 892         }
 893       }
 894     }
 895   }
 896
 897   return aResult;
 898 }
 899
 900 //=============================================================================
 901 /*!
 902  *  Generate the internal group of a Pipe T-Shape
 903  */
 904 //=============================================================================
 905 bool GEOMImpl_IAdvancedOperations::MakeInternalGroup
 906                       (const Handle(GEOM_Object) &theShape,
 907                        const double theR1, const double theLen1,
 908                        const double theR2, const double theLen2,
 909                        const double theRL, double theTransLenL,
 910                        const double theRR, double theTransLenR,
 911                        const double theRI, double theTransLenI,
 912                        const Handle(TColStd_HSequenceOfTransient) &theSeq,
 913                        const gp_Trsf &theTrsf)
 914 {
 915   SetErrorCode(KO);
 916
 917   if (theShape.IsNull()) {
 918     return false;
 919   }
 920
 921   TopoDS_Shape aShape = theShape->GetValue();
 922
 923   if (aShape.IsNull()) {
 924     SetErrorCode("Shape is not defined");
 925     return false;
 926   }
 927
 928   // Compute tolerance
 929   Standard_Real aMaxTol = -RealLast();
 930   TopExp_Explorer anExp(aShape, TopAbs_VERTEX);
 931
 932   for (; anExp.More(); anExp.Next()) {
 933     TopoDS_Vertex aVertex = TopoDS::Vertex(anExp.Current());
 934
 935     if (aVertex.IsNull() == Standard_False) {
 936       const Standard_Real aTol = BRep_Tool::Tolerance(aVertex);
 937
 938       if (aTol > aMaxTol) {
 939         aMaxTol = aTol;
 940       }
 941     }
 942   }
 943
 944   // Construct internal surfaces.
 945   Standard_Integer i = 0;
 946   const Standard_Integer aMaxNbSurf = 5;
 947   Handle(Geom_Surface) aSurface[aMaxNbSurf];
 948   TopTools_ListOfShape aConicalFaces;
 949   Standard_Real aTolConf = Precision::Confusion();
 950
 951   // 1. Construct the internal surface of main pipe.
 952   gp_Ax2 anAxis1 (gp::Origin(), gp::DX(), gp::DZ());
 953   gp_Ax2 anAxis2 (gp::Origin(), gp::DZ(), gp::DX());
 954
 955   aSurface[i++] = new Geom_CylindricalSurface(anAxis1, theR1);
 956
 957   // 2. Construct the internal surface of incident pipe.
 958   aSurface[i++] = new Geom_CylindricalSurface(anAxis2, theR2);
 959
 960   // 3. Construct the internal surface of left reduction pipe.
 961   if (theRL > aTolConf) {
 962     aSurface[i++] = new Geom_CylindricalSurface(anAxis1, theRL);
 963
 964     if (theTransLenL > aTolConf) {
 965       // 3.1. Construct the internal surface of left transition pipe.
 966       gp_Pnt aPLeft (-theLen1, 0., 0.);
 967       gp_Ax2 anAxisLeft (aPLeft, -gp::DX(), gp::DZ());
 968       TopoDS_Shape aConeLeft =
 969         MakeConicalFace(anAxisLeft, theR1, theRL, theTransLenL, theTrsf);
 970
 971       if (aConeLeft.IsNull() == Standard_False) {
 972         aConicalFaces.Append(aConeLeft);
 973       }
 974     }
 975   }
 976
 977   // 4. Construct the internal surface of right reduction pipe.
 978   if (theRR > aTolConf) {
 979     // There is no need to construct another cylinder of the same radius. Skip it.
 980     if (Abs(theRR - theRL) > aTolConf) {
 981       aSurface[i++] = new Geom_CylindricalSurface(anAxis1, theRR);
 982     }
 983
 984     if (theTransLenL > aTolConf) {
 985       // 4.1. Construct the internal surface of right transition pipe.
 986       gp_Pnt aPRight (theLen1, 0., 0.);
 987       gp_Ax2 anAxisRight (aPRight, gp::DX(), gp::DZ());
 988       TopoDS_Shape aConeRight =
 989         MakeConicalFace(anAxisRight, theR1, theRR, theTransLenR, theTrsf);
 990
 991       if (aConeRight.IsNull() == Standard_False) {
 992         aConicalFaces.Append(aConeRight);
 993       }
 994     }
 995   }
 996
 997   // 5. Construct the internal surface of incident reduction pipe.
 998   if (theRI > aTolConf) {
 999     aSurface[i++] = new Geom_CylindricalSurface(anAxis2, theRI);
1000
1001     if (theTransLenI > aTolConf) {
1002       // 5.1. Construct the internal surface of incident transition pipe.
1003       gp_Pnt aPInci (0., 0., theLen2);
1004       gp_Ax2 anAxisInci (aPInci, gp::DZ(), gp::DX());
1005       TopoDS_Shape aConeInci =
1006         MakeConicalFace(anAxisInci, theR2, theRI, theTransLenI, theTrsf);
1007
1008       if (aConeInci.IsNull() == Standard_False) {
1009         aConicalFaces.Append(aConeInci);
1010       }
1011     }
1012   }
1013
1014   // Get faces that are laying on cylindrical surfaces.
1015   TopTools_ListOfShape aFaces;
1016   gp_Trsf anInvTrsf = theTrsf.Inverted();
1017
1018   for (i = 0; i < aMaxNbSurf; i++) {
1019     if (aSurface[i].IsNull()) {
1020       break;
1021     }
1022
1023     aSurface[i]->Transform(theTrsf);
1024
1025     TopTools_ListOfShape aLocalFaces;
1026
1027     if (!GetFacesOnSurf(aShape, aSurface[i], aMaxTol, aLocalFaces)) {
1028       return false;
1029     }
1030
1031     if (i < 2) {
1032       // Check if the result contains outer cylinders.
1033       // It is required for main and incident pipes.
1034       TopTools_ListIteratorOfListOfShape anIter(aLocalFaces);
1035
1036       while (anIter.More()) {
1037         TopExp_Explorer anExp(anIter.Value(), TopAbs_VERTEX);
1038         Standard_Boolean isInside = Standard_False;
1039
1040         // Get a vertex from this shape
1041         if (anExp.More()) {
1042           TopoDS_Vertex aVtx = TopoDS::Vertex(anExp.Current());
1043
1044           if (aVtx.IsNull() == Standard_False) {
1045             gp_Pnt aPnt = BRep_Tool::Pnt(aVtx);
1046
1047             aPnt.Transform(anInvTrsf);
1048
1049             if (i == 0) {
1050               // Check if the point is inside the main pipe.
1051               isInside = (Abs(aPnt.X()) <= theLen1);
1052             } else { // i == 1
1053               // Check if the point is inside the incident pipe.
1054               isInside = (aPnt.Z() <= theLen2);
1055             }
1056           }
1057         }
1058
1059         if (isInside) {
1060           // Keep this face.
1061           anIter.Next();
1062         } else {
1063           // Remove this face.
1064           aLocalFaces.Remove(anIter);
1065         }
1066       }
1067     }
1068
1069     aFaces.Append(aLocalFaces);
1070   }
1071
1072   // Get faces that are laying on conical faces.
1073   if (aConicalFaces.IsEmpty() == Standard_False) {
1074     Handle(GEOM_Object) aCone =
1075       GetEngine()->AddObject(GetDocID(), GEOM_TSHAPE);
1076     Handle(GEOM_Function) aFunction =
1077       aCone->AddFunction(GEOMImpl_PipeTShapeDriver::GetID(), TSHAPE_BASIC);
1078     TopTools_ListIteratorOfListOfShape aFIter(aConicalFaces);
1079     Handle(GEOM_Object) aConeFromShape;
1080
1081     for (; aFIter.More(); aFIter.Next()) {
1082       aFunction->SetValue(aFIter.Value());
1083       aConeFromShape = myShapesOperations->GetInPlace(theShape, aCone);
1084
1085       if (aConeFromShape.IsNull() == Standard_False) {
1086         aConeFromShape->GetLastFunction()->SetDescription("");
1087         TopoDS_Shape aConeFaces = aConeFromShape->GetValue();
1088         TopExp_Explorer anExp(aConeFaces, TopAbs_FACE);
1089
1090         for (; anExp.More(); anExp.Next()) {
1091           TopoDS_Face aConeFace = TopoDS::Face(anExp.Current());
1092
1093           if (aConeFace.IsNull() == Standard_False) {
1094             aFaces.Append(aConeFace);
1095           }
1096         }
1097       }
1098     }
1099   }
1100
1101   // Create a group of internal faces.
1102   if (aFaces.IsEmpty() == Standard_False) {
1103     Handle(GEOM_Object) aGroup = myGroupOperations->CreateGroup(theShape, TopAbs_FACE);
1104
1105     if (aGroup.IsNull() == Standard_False) {
1106       aGroup->GetLastFunction()->SetDescription("");
1107       aGroup->SetName("INTERNAL_FACES");
1108
1109       TopTools_IndexedMapOfShape anIndices;
1110       Handle(TColStd_HSequenceOfInteger) aSeqIDs = new TColStd_HSequenceOfInteger;
1111
1112       TopExp::MapShapes(aShape, anIndices);
1113
1114       TopTools_ListIteratorOfListOfShape anIter(aFaces);
1115
1116       for (; anIter.More(); anIter.Next()) {
1117         const TopoDS_Shape &aFace = anIter.Value();
1118         const Standard_Integer anIndex = anIndices.FindIndex(aFace);
1119
1120         if (anIndex > 0) {
1121           aSeqIDs->Append(anIndex);
1122         }
1123       }
1124
1125       myGroupOperations->UnionIDs(aGroup, aSeqIDs);
1126       aGroup->GetLastFunction()->SetDescription("");
1127       theSeq->Append(aGroup);
1128     }
1129   }
1130
1131   SetErrorCode(OK);
1132
1133   return true;
1134 }
1135
1136 bool GEOMImpl_IAdvancedOperations::MakePipeTShapePartition(Handle(GEOM_Object) theShape,
1137                                                            double theR1, double theW1, double theL1,
1138                                                            double theR2, double theW2, double theL2,
1139                                                            double theH, double theW,
1140                                                            double theRF, bool isNormal)
1141 {
1142   SetErrorCode(KO);
1143
1144   // Build tools for partition operation:
1145   // 1 face and 2 planes
1146   // Face
1147   Handle(GEOM_Object) arete_intersect_int, arete_intersect_ext;
1148   Handle(GEOM_Object) wire_t, wire_t2, face_t, face_t2;
1149   Handle(GEOM_Object) chan_racc;
1150   Handle(GEOM_Object) vi1, vi2;
1151   Handle(GEOM_Object) Te3;
1152
1153   try {
1154 #if OCC_VERSION_LARGE > 0x06010000
1155     OCC_CATCH_SIGNALS;
1156 #endif
1157     Handle(GEOM_Object) Vector_Z = myBasicOperations->MakeVectorDXDYDZ(0, 0, 1);
1158     Vector_Z->GetLastFunction()->SetDescription("");
1159
1160     // Useful values
1161     double aSize = 2*(theL1 + theL2);
1162     double aR1Ext = theR1 + theW1;
1163     double aR2Ext = theR2 + theW2;
1164     double theVertCylinderRadius = aR2Ext + theW + theRF;
1165     double theHoriCylinderRadius = aR1Ext + theH + theRF;
1166
1167     // Common edges on internal cylinder
1168     Handle(GEOM_Object) box_i = my3DPrimOperations->MakeBoxDXDYDZ(theR2, theR2, theR1);
1169     box_i->GetLastFunction()->SetDescription("");
1170     box_i = myTransformOperations->TranslateDXDYDZ(box_i, -theR2, -theR2, 0);
1171     box_i->GetLastFunction()->SetDescription("");
1172
1173     Handle(GEOM_Function) aFunction = theShape->GetLastFunction();
1174     TCollection_AsciiString theDesc = aFunction->GetDescription();
1175     Handle(TColStd_HSequenceOfTransient) edges_i =
1176       myShapesOperations->GetShapesOnBox(box_i, theShape, TopAbs_EDGE, GEOMAlgo_ST_IN);
1177     // Recover previous description to get rid of Propagate dump
1178     aFunction->SetDescription(theDesc);
1179     if (edges_i.IsNull() || edges_i->Length() == 0) {
1180       SetErrorCode("Internal edges not found");
1181       return false;
1182     }
1183     for (int i=1; i<=edges_i->Length();i++) {
1184       Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast(edges_i->Value(i));
1185       anObj->GetLastFunction()->SetDescription("");
1186     }
1187     arete_intersect_int = Handle(GEOM_Object)::DownCast(edges_i->Value(1));
1188
1189     // search for vertices located on both internal pipes
1190     aFunction = theShape->GetLastFunction();
1191     theDesc = aFunction->GetDescription();
1192     Handle(TColStd_HSequenceOfTransient) vertices_i =
1193       myShapesOperations->GetShapesOnBox(box_i, theShape, TopAbs_VERTEX, GEOMAlgo_ST_ONIN);
1194     // Recover previous description to get rid of Propagate dump
1195     aFunction->SetDescription(theDesc);
1196     if (vertices_i.IsNull() || vertices_i->Length() == 0) {
1197       SetErrorCode("Internal vertices not found");
1198       return false;
1199     }
1200
1201     double d1min = theR2+theW2, d2min=theR2+theW2;
1202     for (int i = 1; i <= vertices_i->Length(); i++) {
1203       Handle(GEOM_Object) v = Handle(GEOM_Object)::DownCast(vertices_i->Value(i));
1204       v->GetLastFunction()->SetDescription("");
1205       TopoDS_Vertex aVertex = TopoDS::Vertex(v->GetValue());
1206       gp_Pnt aP = BRep_Tool::Pnt(aVertex);
1207       if (Abs(aP.X()) <= Precision::Confusion()) {
1208         if (Abs(aP.Y()) < d1min) {
1209           vi1 = v;
1210           d1min = Abs(aP.Y());
1211         }
1212       } else if (Abs(aP.Y()) <= Precision::Confusion()) {
1213         if (Abs(aP.X()) < d2min) {
1214           vi2 = v;
1215           d2min = Abs(aP.X());
1216         }
1217       }
1218     }
1219     if (vi1.IsNull() || vi2.IsNull()) {
1220       SetErrorCode("Cannot find internal intersection vertices");
1221       return false;
1222     }
1223
1224     std::list<Handle(GEOM_Object)> theShapes;
1225
1226     if (isNormal) {
1227       Handle(GEOM_Object) ve1, ve2;
1228       TopoDS_Vertex vertex1, vertex2;
1229
1230       Handle(GEOM_Object) box_e = my3DPrimOperations->MakeBoxDXDYDZ(aR2Ext, aR2Ext, aR1Ext);
1231       box_e->GetLastFunction()->SetDescription("");
1232       box_e = myTransformOperations->TranslateDXDYDZ(box_e, -aR2Ext, -aR2Ext, 0);
1233       box_e->GetLastFunction()->SetDescription("");
1234
1235       // search for vertices located on both external pipes
1236       aFunction = theShape->GetLastFunction();
1237       theDesc = aFunction->GetDescription();
1238       Handle(TColStd_HSequenceOfTransient) vertices_e =
1239         myShapesOperations->GetShapesOnBox(box_e, theShape, TopAbs_VERTEX, GEOMAlgo_ST_ONIN);
1240       // Recover previous description to get rid of Propagate dump
1241       aFunction->SetDescription(theDesc);
1242       if (vertices_e.IsNull() || vertices_e->Length() == 0) {
1243         SetErrorCode("External vertices not found");
1244         return false;
1245       }
1246
1247       double d1max = 0, d2max = 0;
1248       for (int i = 1; i <= vertices_e->Length(); i++) {
1249         Handle(GEOM_Object) v = Handle(GEOM_Object)::DownCast(vertices_e->Value(i));
1250         v->GetLastFunction()->SetDescription("");
1251         TopoDS_Vertex aVertex = TopoDS::Vertex(v->GetValue());
1252         gp_Pnt aP = BRep_Tool::Pnt(aVertex);
1253         if (Abs(aP.X()) <= Precision::Confusion()) {
1254           if (Abs(aP.Y()) > d1max) {
1255             ve1 = v;
1256             vertex1 = aVertex;
1257             d1max = Abs(aP.Y());
1258           }
1259         } else if (Abs(aP.Y()) <= Precision::Confusion()) {
1260           if (Abs(aP.X()) > d2max) {
1261             ve2 = v;
1262             vertex2 = aVertex;
1263             d2max = Abs(aP.X());
1264           }
1265         }
1266       }
1267       if (ve1.IsNull() || ve2.IsNull()) {
1268         SetErrorCode("Cannot find external intersection vertices");
1269         return false;
1270       }
1271       Handle(GEOM_Object) edge_e1, edge_e2;
1272
1273       // Common edges on external cylinder
1274       aFunction = theShape->GetLastFunction();
1275       theDesc = aFunction->GetDescription();
1276       Handle(TColStd_HSequenceOfTransient) edges_e =
1277         myShapesOperations->GetShapesOnBox(box_e, theShape, TopAbs_EDGE, GEOMAlgo_ST_IN);
1278       // Recover previous description to get rid of Propagate dump
1279       aFunction->SetDescription(theDesc);
1280       if (edges_e.IsNull() || edges_e->Length() == 0) {
1281         SetErrorCode("External edges not found");
1282         return false;
1283       }
1284       ShapeAnalysis_Edge sae;
1285       for (int i=1; i<=edges_e->Length();i++) {
1286         Handle(GEOM_Object) anObj = Handle(GEOM_Object)::DownCast(edges_e->Value(i));
1287         anObj->GetLastFunction()->SetDescription("");
1288         TopoDS_Edge anEdge = TopoDS::Edge(anObj->GetValue());
1289         if ( !anEdge.IsNull() &&
1290              (sae.FirstVertex(anEdge).IsSame(vertex1) || sae.LastVertex(anEdge).IsSame(vertex1)) &&
1291              (sae.FirstVertex(anEdge).IsSame(vertex2) || sae.LastVertex(anEdge).IsSame(vertex2))) {
1292           arete_intersect_ext = anObj;
1293         }
1294       }
1295
1296       edge_e1 = myBasicOperations->MakeLineTwoPnt(ve1, vi1);
1297       if (edge_e1.IsNull()) {
1298         SetErrorCode("Edge 1 could not be built");
1299         return false;
1300       }
1301
1302       edge_e2 = myBasicOperations->MakeLineTwoPnt(ve2, vi2);
1303       if (edge_e2.IsNull()) {
1304         SetErrorCode("Edge 2 could not be built");
1305         return false;
1306       }
1307
1308       edge_e1->GetLastFunction()->SetDescription("");
1309       edge_e2->GetLastFunction()->SetDescription("");
1310
1311       std::list<Handle(GEOM_Object)> edge_e_elist;
1312       edge_e_elist.push_back(arete_intersect_int);
1313       edge_e_elist.push_back(edge_e1);
1314       edge_e_elist.push_back(arete_intersect_ext);
1315       edge_e_elist.push_back(edge_e2);
1316       wire_t = myShapesOperations->MakeWire(edge_e_elist, 1e-7);
1317       if (wire_t.IsNull()) {
1318         SetErrorCode("Impossible to build wire");
1319         return false;
1320       }
1321       wire_t->GetLastFunction()->SetDescription("");
1322       face_t = myShapesOperations->MakeFace(wire_t, false);
1323       if (face_t.IsNull()) {
1324         SetErrorCode("Impossible to build face");
1325         return false;
1326       }
1327       face_t->GetLastFunction()->SetDescription("");
1328
1329       theShapes.push_back(theShape);
1330       theShapes.push_back(vi1);
1331       theShapes.push_back(vi2);
1332       theShapes.push_back(ve1);
1333       theShapes.push_back(ve2);
1334       theShapes.push_back(edge_e1);
1335       theShapes.push_back(edge_e2);
1336       theShapes.push_back(wire_t);
1337       theShapes.push_back(face_t);
1338     }
1339     else {
1340       Handle(GEOM_Object) P1, P2, P3, P4, P5, P6;
1341       int idP1, idP2, idP3, idP4;
1342       int PZX, PZY;
1343       double ZX=0, ZY=0;
1344       std::vector<int> LX;
1345       std::vector<int> LY;
1346       Handle(GEOM_Object) box_e = my3DPrimOperations->MakeBoxDXDYDZ
1347         (theVertCylinderRadius, theVertCylinderRadius, theHoriCylinderRadius);
1348       box_e->GetLastFunction()->SetDescription("");
1349       box_e = myTransformOperations->TranslateDXDYDZ
1350         (box_e, -theVertCylinderRadius, -theVertCylinderRadius, 0);
1351       box_e->GetLastFunction()->SetDescription("");
1352
1353       aFunction = theShape->GetLastFunction();
1354       theDesc = aFunction->GetDescription();
1355       Handle(TColStd_HSequenceOfTransient) extremVertices =
1356         myShapesOperations->GetShapesOnBox(box_e, theShape, TopAbs_VERTEX, GEOMAlgo_ST_ONIN);
1357       // Recover previous description to get rid of Propagate dump
1358       aFunction->SetDescription(theDesc);
1359
1360       if (extremVertices.IsNull() || extremVertices->Length() == 0) {
1361         if (theRF == 0)
1362           SetErrorCode("Vertices on chamfer not found");
1363         else
1364           SetErrorCode("Vertices on fillet not found");
1365         return false;
1366       }
1367
1368       theShapes.push_back(theShape);
1369       theShapes.push_back(box_e);
1370       if (extremVertices->Length() != 6) {
1371         //           for (int i=1; i<=extremVertices->Length(); i++){
1372         //             theShapes.push_back(Handle(GEOM_Object)::DownCast(extremVertices->Value(i)));
1373         //           }
1374         //           Handle(GEOM_Object) aCompound = myShapesOperations->MakeCompound(theShapes);
1375         //           TopoDS_Shape aCompoundShape = aCompound->GetValue();
1376         //           theShape->GetLastFunction()->SetValue(aCompoundShape);
1377         SetErrorCode("Bad number of vertices on chamfer found");
1378         return false;
1379       }
1380
1381       for (int i=1; i<=extremVertices->Length(); i++){
1382         Handle(GEOM_Object) aV = Handle(GEOM_Object)::DownCast(extremVertices->Value(i));
1383         aV->GetLastFunction()->SetDescription("");
1384         gp_Pnt aP = BRep_Tool::Pnt(TopoDS::Vertex(aV->GetValue()));
1385
1386         if (Abs(aP.X()) <= Precision::Confusion()) {
1387           if (Abs(aP.Y()) - theR2 > Precision::Confusion()) {
1388             LX.push_back(i);
1389             if  (aP.Z()-ZX > Precision::Confusion()) {
1390               ZX = aP.Z();
1391               PZX = i;
1392             }
1393           }
1394         }
1395         else {
1396           if (Abs(aP.X()) - theR2 > Precision::Confusion()) {
1397             LY.push_back(i);
1398             if (aP.Z() - ZY > Precision::Confusion()) {
1399               ZY = aP.Z();
1400               PZY = i;
1401             }
1402           }
1403         }
1404       }
1405
1406       idP2 = PZX;
1407       idP4 = PZY;
1408       idP1 = LX.at(0);
1409       if (LX.at(0) == PZX)
1410         idP1 = LX.at(1);
1411       idP3 = LY.at(0);
1412       if (LY.at(0) == PZY)
1413         idP3 = LY.at(1);
1414
1415       P1 = Handle(GEOM_Object)::DownCast(extremVertices->Value(idP1));
1416       P2 = Handle(GEOM_Object)::DownCast(extremVertices->Value(idP2));
1417       P3 = Handle(GEOM_Object)::DownCast(extremVertices->Value(idP3));
1418       P4 = Handle(GEOM_Object)::DownCast(extremVertices->Value(idP4));
1419
1420       Handle(GEOM_Object) Cote_1 = myBasicOperations->MakeLineTwoPnt(P1, vi1);
1421       if (Cote_1.IsNull()) {
1422         SetErrorCode("Impossible to build edge in thickness");
1423         return false;
1424       }
1425       Cote_1->GetLastFunction()->SetDescription("");
1426
1427       Handle(GEOM_Object) Cote_2 = myBasicOperations->MakeLineTwoPnt(vi2, P3);
1428       if (Cote_2.IsNull()) {
1429         SetErrorCode("Impossible to build edge in thickness");
1430         return false;
1431       }
1432       Cote_2->GetLastFunction()->SetDescription("");
1433
1434       // edge_chan_princ = arete du chanfrein (ou raccord) sur le tuyau principal
1435       // edge_chan_inc = arete du chanfrein (ou raccord) sur le tuyau incident
1436       //         std::cerr << "Getting chamfer edge on main pipe" << std::endl;
1437       Handle(GEOM_Object) edge_chan_princ = myBlocksOperations->GetEdge(theShape, P1, P3);
1438       if (edge_chan_princ.IsNull()) {
1439         SetErrorCode("Impossible to find edge on main pipe");
1440         return false;
1441       }
1442       edge_chan_princ->GetLastFunction()->SetDescription("");
1443
1444       Handle(GEOM_Object) edge_chan_inc = myBlocksOperations->GetEdge(theShape, P2, P4);
1445       if (edge_chan_inc.IsNull()) {
1446         SetErrorCode("Impossible to find edge on incident pipe");
1447         return false;
1448       }
1449       edge_chan_inc->GetLastFunction()->SetDescription("");
1450
1451       std::list<Handle(GEOM_Object)> edgeList1;
1452       edgeList1.push_back(edge_chan_princ);
1453       edgeList1.push_back(Cote_1);
1454       edgeList1.push_back(arete_intersect_int);
1455       edgeList1.push_back(Cote_2);
1456
1457       //         std::cerr << "Creating wire 1" << std::endl;
1458       wire_t = myShapesOperations->MakeWire(edgeList1, 1e-7);
1459       if (wire_t.IsNull()) {
1460         SetErrorCode("Impossible to build wire");
1461         return false;
1462       }
1463       wire_t->GetLastFunction()->SetDescription("");
1464
1465       //         std::cerr << "Creating face 1" << std::endl;
1466       face_t = myShapesOperations->MakeFace(wire_t, false);
1467       if (face_t.IsNull()) {
1468         SetErrorCode("Impossible to build face");
1469         return false;
1470       }
1471       face_t->GetLastFunction()->SetDescription("");
1472       theShapes.push_back(face_t);
1473
1474       gp_Pnt aP2 = BRep_Tool::Pnt(TopoDS::Vertex(P2->GetValue()));
1475       gp_Pnt aP5 = BRep_Tool::Pnt(TopoDS::Vertex(vi1->GetValue()));
1476       double deltaZ = aP2.Z() - aP5.Z();
1477       //         std::cerr << "Creating new point from vi1 with deltaZ = " << deltaZ << std::endl;
1478       Handle(GEOM_Object) P5bis = myTransformOperations->TranslateDXDYDZCopy(vi1, 0, 0, deltaZ);
1479       if (P5bis.IsNull()) {
1480         SetErrorCode("Impossible to translate vertex");
1481         return false;
1482       }
1483       P5bis->GetLastFunction()->SetDescription("");
1484
1485       gp_Pnt aP4 = BRep_Tool::Pnt(TopoDS::Vertex(P4->GetValue()));
1486       gp_Pnt aP6 = BRep_Tool::Pnt(TopoDS::Vertex(vi2->GetValue()));
1487       deltaZ = aP4.Z() - aP6.Z();
1488       //         std::cerr << "Creating new point from vi2 with deltaZ = " << deltaZ << std::endl;
1489       Handle(GEOM_Object) P6bis = myTransformOperations->TranslateDXDYDZCopy(vi2, 0, 0, deltaZ);
1490       if (P6bis.IsNull()) {
1491         SetErrorCode("Impossible to translate vertex");
1492         return false;
1493       }
1494       P6bis->GetLastFunction()->SetDescription("");
1495
1496       //         std::cerr << "Creating new line 1 from 2 previous points" << std::endl;
1497       Handle(GEOM_Object) Cote_3 = myBasicOperations->MakeLineTwoPnt(P5bis, P2);
1498       if (Cote_3.IsNull()) {
1499         SetErrorCode("Impossible to build edge in thickness");
1500         return false;
1501       }
1502       Cote_3->GetLastFunction()->SetDescription("");
1503
1504       //         std::cerr << "Creating new line 2 from 2 previous points" << std::endl;
1505       Handle(GEOM_Object) Cote_4 = myBasicOperations->MakeLineTwoPnt(P6bis, P4);
1506       if (Cote_4.IsNull()) {
1507         SetErrorCode("Impossible to build edge in thickness");
1508         return false;
1509       }
1510       Cote_4->GetLastFunction()->SetDescription("");
1511
1512       //         std::cerr << "Creating new line 3 from 2 previous points" << std::endl;
1513       Handle(GEOM_Object) Cote_5 = myBasicOperations->MakeLineTwoPnt(P5bis, P6bis);
1514       if (Cote_4.IsNull()) {
1515         SetErrorCode("Impossible to build edge in thickness");
1516         return false;
1517       }
1518       Cote_5->GetLastFunction()->SetDescription("");
1519
1520       //std::list<Handle(GEOM_Object)> edgeList2;
1521       //edgeList2.push_back(edge_chan_inc);
1522       //edgeList2.push_back(Cote_3);
1523       //edgeList2.push_back(Cote_5);
1524       //edgeList2.push_back(Cote_4);
1525       //         std::cerr << "Creating wire 2" << std::endl;
1526       //wire_t2 = myShapesOperations->MakeWire(edgeList2, 1e-7);
1527       //if (wire_t2.IsNull()) {
1528       //  SetErrorCode("Impossible to build wire");
1529       //  return false;
1530       //}
1531       //wire_t2->GetLastFunction()->SetDescription("");
1532       //         std::cerr << "Creating face 2" << std::endl;
1533       //face_t2 = myShapesOperations->MakeFace(wire_t2, false);
1534
1535       // Mantis issue 0021682
1536       face_t2 = my3DPrimOperations->MakePrismVecH(edge_chan_inc, Cote_4, - (theR2 + theW2));
1537       //face_t2 = my3DPrimOperations->MakePrismVecH(edge_chan_inc, Cote_4, - 2.0*theR2);
1538       if (face_t2.IsNull()) {
1539         SetErrorCode("Impossible to build face");
1540         return false;
1541       }
1542       face_t2->GetLastFunction()->SetDescription("");
1543       theShapes.push_back(face_t2);
1544     }
1545
1546     // Planes
1547     Handle(GEOM_Object) aP0 = myBasicOperations->MakePointXYZ(0, 0, 0);
1548     Handle(GEOM_Object) aVZ = myBasicOperations->MakeVectorDXDYDZ(0, 0, 1);
1549     Handle(GEOM_Object) aVXZ = myBasicOperations->MakeVectorDXDYDZ(aR1Ext, 0, 0.5*(theL1+theVertCylinderRadius));
1550     Handle(GEOM_Object) aPlnOZ = myBasicOperations->MakePlanePntVec(aP0, aVZ, aSize);
1551     Handle(GEOM_Object) aPlnOXZ = myBasicOperations->MakePlanePntVec(aP0, aVXZ, aSize);
1552     aP0->GetLastFunction()->SetDescription("");
1553     aVZ->GetLastFunction()->SetDescription("");
1554     aVXZ->GetLastFunction()->SetDescription("");
1555     aPlnOZ->GetLastFunction()->SetDescription("");
1556     aPlnOXZ->GetLastFunction()->SetDescription("");
1557     theShapes.push_back(aPlnOZ);
1558     theShapes.push_back(aPlnOXZ);
1559
1560     // Partition
1561     Handle(TColStd_HSequenceOfTransient) partitionShapes = new TColStd_HSequenceOfTransient;
1562     Handle(TColStd_HSequenceOfTransient) theTools = new TColStd_HSequenceOfTransient;
1563     Handle(TColStd_HSequenceOfTransient) theKeepInside = new TColStd_HSequenceOfTransient;
1564     Handle(TColStd_HSequenceOfTransient) theRemoveInside = new TColStd_HSequenceOfTransient;
1565     Handle(TColStd_HArray1OfInteger) theMaterials;
1566
1567     partitionShapes->Append(theShape);
1568     theTools->Append(aPlnOZ);
1569     if (Abs(aR1Ext - aR2Ext) > Precision::Confusion())
1570       theTools->Append(aPlnOXZ);
1571     theTools->Append(face_t);
1572     if (!isNormal)
1573       theTools->Append(face_t2);
1574
1575     Te3 = myBooleanOperations->MakePartition
1576               (partitionShapes, theTools, theKeepInside, theRemoveInside,
1577               TopAbs_SOLID, false, theMaterials, 0, false);
1578     if (Te3.IsNull()) {
1579       SetErrorCode("Impossible to build partition of TShape");
1580       return false;
1581     }
1582     Te3->GetLastFunction()->SetDescription("");
1583
1584     // Last verification: result should be a block
1585     std::list<GEOMImpl_IBlocksOperations::BCError> errList;
1586     if (!myBlocksOperations->CheckCompoundOfBlocks(Te3,errList)) {
1587       SetErrorCode("TShape is not a compound of block");
1588       return false;
1589     }
1590
1591 //     // BEGIN Compound of created shapes - Only for debug purpose
1592 //     theShapes.clear();
1593 //     theShapes.push_back(theShape);
1594 //     theShapes.push_back(aPlnOZ);
1595 //     if (Abs(aR1Ext - aR2Ext) > Precision::Confusion() )
1596 //       theShapes.push_back(aPlnOXZ);
1597 //     theShapes.push_back(face_t);
1598 //     if (!isNormal)
1599 //       theShapes.push_back(face_t2);
1600 //
1601 //     Handle(GEOM_Object) aCompound = myShapesOperations->MakeCompound(theShapes);
1602 //     TopoDS_Shape aCompoundShape = aCompound->GetValue();
1603 //     theShape->GetLastFunction()->SetValue(aCompoundShape);
1604 //     // END Compound of created shapes - Only for debug purpose
1605
1606     TopoDS_Shape aShape = Te3->GetValue();
1607     theShape->GetLastFunction()->SetValue(aShape);
1608   }
1609   catch (Standard_Failure) {
1610     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
1611     SetErrorCode(aFail->GetMessageString());
1612     return false;
1613   }
1614
1615   SetErrorCode(OK);
1616   return true;
1617 }
1618
1619 // Mirror and glue faces
1620 bool GEOMImpl_IAdvancedOperations::MakePipeTShapeMirrorAndGlue(Handle(GEOM_Object) theShape,
1621                                                                double theR1, double theW1, double theL1,
1622                                                                double theR2, double theW2, double theL2)
1623 {
1624   SetErrorCode(KO);
1625
1626   // Useful values
1627   double aSize = 2*(theL1 + theL2);
1628   double aR1Ext = theR1 + theW1;
1629
1630   // Planes
1631   Handle(GEOM_Object) aP0 = myBasicOperations->MakePointXYZ(0, 0, 0);
1632   aP0->GetLastFunction()->SetDescription("");
1633   Handle(GEOM_Object) aVX = myBasicOperations->MakeVectorDXDYDZ(1, 0, 0);
1634   Handle(GEOM_Object) aVY = myBasicOperations->MakeVectorDXDYDZ(0, 1, 0);
1635   aVX->GetLastFunction()->SetDescription("");
1636   aVY->GetLastFunction()->SetDescription("");
1637   Handle(GEOM_Object) aPlane_OX = myBasicOperations->MakePlanePntVec(aP0, aVX, 2*(aR1Ext + theL2));
1638   Handle(GEOM_Object) aPlane_OY = myBasicOperations->MakePlanePntVec(aP0, aVY, aSize);
1639   aPlane_OX->GetLastFunction()->SetDescription("");
1640   aPlane_OY->GetLastFunction()->SetDescription("");
1641
1642   Handle(GEOM_Object) Te4 = myTransformOperations->MirrorPlaneCopy(theShape, aPlane_OX);
1643   if (Te4.IsNull()) {
1644     SetErrorCode("Impossible to build mirror of quarter TShape");
1645     return false;
1646   }
1647
1648   Handle(GEOM_Object) Te5 = myTransformOperations->MirrorPlaneCopy(theShape, aPlane_OY);
1649   if (Te5.IsNull()) {
1650     SetErrorCode("Impossible to build mirror of half TShape");
1651     return false;
1652   }
1653
1654   Handle(GEOM_Object) Te6 = myTransformOperations->MirrorPlaneCopy(Te4, aPlane_OY);
1655   if (Te6.IsNull()) {
1656     SetErrorCode("Impossible to build mirror of half TShape");
1657     return false;
1658   }
1659
1660   std::list<Handle(GEOM_Object)> aShapesList;
1661   aShapesList.push_back(theShape);
1662   aShapesList.push_back(Te4);
1663   aShapesList.push_back(Te5);
1664   aShapesList.push_back(Te6);
1665   Handle(GEOM_Object) Te7 = myShapesOperations->MakeCompound(aShapesList);
1666   if (Te7.IsNull()) {
1667     SetErrorCode("Impossible to build compound");
1668     return false;
1669   }
1670
1671   // Copy source shape
1672   TopoDS_Shape aShapeCopy;
1673   TColStd_IndexedDataMapOfTransientTransient aMapTShapes;
1674   TNaming_CopyShape::CopyTool(Te7->GetValue(), aMapTShapes, aShapeCopy);
1675
1676   Handle(GEOM_Object) Te8 = myShapesOperations->MakeGlueFaces(Te7, 1e-7, true);
1677   if (Te8.IsNull()) {
1678     SetErrorCode("Impossible to glue faces of TShape");
1679     return false;
1680   }
1681
1682   TopoDS_Shape aShape = Te8->GetValue();
1683   BRepCheck_Analyzer anAna (aShape, Standard_True);
1684
1685   if (!anAna.IsValid()) {
1686     // Try to do gluing with the tolerance equal to maximal
1687     // tolerance of vertices of the source shape.
1688     Standard_Real aTolMax = -RealLast();
1689
1690     for (TopExp_Explorer ExV (aShapeCopy, TopAbs_VERTEX); ExV.More(); ExV.Next()) {
1691       TopoDS_Vertex aVertex = TopoDS::Vertex(ExV.Current());
1692       Standard_Real aTol = BRep_Tool::Tolerance(aVertex);
1693
1694       if (aTol > aTolMax) {
1695         aTolMax = aTol;
1696       }
1697     }
1698
1699     // Perform gluing
1700     Te7->GetLastFunction()->SetValue(aShapeCopy);
1701     Te8 = myShapesOperations->MakeGlueFaces(Te7, aTolMax, true);
1702
1703     if (Te8.IsNull()) {
1704       SetErrorCode("Impossible to glue faces of TShape");
1705       return false;
1706     }
1707
1708     aShape = Te8->GetValue();
1709   }
1710
1711
1712   theShape->GetLastFunction()->SetValue(aShape);
1713
1714   Te4->GetLastFunction()->SetDescription("");
1715   Te5->GetLastFunction()->SetDescription("");
1716   Te6->GetLastFunction()->SetDescription("");
1717   Te7->GetLastFunction()->SetDescription("");
1718   Te8->GetLastFunction()->SetDescription("");
1719
1720   SetErrorCode(OK);
1721   return true;
1722 }
1723
1724 //=======================================================================
1725 //function : MakePipeTShapeThicknessReduction
1726 //purpose  : Static method. Add thiskness reduction elements at the three
1727 //                          open ends of the T-Shape.
1728 //=======================================================================
1729 TopoDS_Shape GEOMImpl_IAdvancedOperations::MakePipeTShapeThicknessReduction
1730                                   (TopoDS_Shape theShape,
1731                                    double r1, double w1, double l1,
1732                                    double r2, double w2, double l2,
1733                                    double rL, double wL, double ltransL, double lthinL,
1734                                    double rR, double wR, double ltransR, double lthinR,
1735                                    double rI, double wI, double ltransI, double lthinI,
1736                                    bool fuseReductions)
1737 {
1738   // Add thickness reduction elements
1739   // at the three extremities: Left, Right and Incident
1740   //
1741   // ---------------------.
1742   //   W                   \
1743   // ---------------------. \
1744   //   ^                   \ '-----------------.
1745   //   |R                   \        Wthin     |
1746   //   |                     '-----------------'
1747   //   v                             Rthin
1748   // --.--.--.--.--.--.--.--.--.--.--.--.--.--.--
1749   //                     Ltrans    Lthin
1750
1751   TopoDS_Shape aResult = theShape;
1752   double aTol = Precision::Confusion();
1753
1754   gp_Vec aVX = gp::DX(), aVZ = gp::DZ();
1755
1756   // Left reduction (rL, wL, ltransL, lthinL)
1757   if (rL > aTol && wL > aTol && ltransL > aTol) {
1758     gp_Pnt aPLeft (-l1, 0, 0);
1759     gp_Ax2 anAxesLeft (aPLeft, -aVX, aVZ);
1760     TopoDS_Shape aReductionLeft = GEOMImpl_IAdvancedOperations::MakeThicknessReduction
1761       (anAxesLeft, r1, w1, rL, wL, ltransL, lthinL, fuseReductions);
1762
1763     if (fuseReductions) {
1764       BRepAlgoAPI_Fuse fuseL (aResult, aReductionLeft);
1765       if (!fuseL.IsDone())
1766         StdFail_NotDone::Raise("Cannot fuse Te with left reduction");
1767       aResult = fuseL.Shape();
1768     }
1769     else {
1770       BRep_Builder B;
1771       TopoDS_Compound C;
1772       B.MakeCompound(C);
1773       B.Add(C, aResult);
1774       B.Add(C, aReductionLeft);
1775       aResult = C;
1776     }
1777   }
1778
1779   // Right reduction
1780   if (rR > aTol && wR > aTol && ltransR > aTol) {
1781     gp_Pnt aPRight (l1, 0, 0);
1782     gp_Ax2 anAxesRight (aPRight, aVX, aVZ);
1783     TopoDS_Shape aReductionRight = GEOMImpl_IAdvancedOperations::MakeThicknessReduction
1784       (anAxesRight, r1, w1, rR, wR, ltransR, lthinR, fuseReductions);
1785
1786     if (fuseReductions) {
1787       BRepAlgoAPI_Fuse fuseR (aResult, aReductionRight);
1788       if (!fuseR.IsDone())
1789         StdFail_NotDone::Raise("Cannot fuse Te with right reduction");
1790       aResult = fuseR.Shape();
1791     }
1792     else {
1793       BRep_Builder B;
1794       TopoDS_Compound C;
1795       B.MakeCompound(C);
1796       B.Add(C, aResult);
1797       B.Add(C, aReductionRight);
1798       aResult = C;
1799     }
1800   }
1801
1802   // Incident reduction
1803   if (rI > aTol && wI > aTol && ltransI > aTol) {
1804     gp_Pnt aPInci (0, 0, l2);
1805     gp_Ax2 anAxesInci (aPInci, aVZ, aVX);
1806     TopoDS_Shape aReductionInci = GEOMImpl_IAdvancedOperations::MakeThicknessReduction
1807       (anAxesInci, r2, w2, rI, wI, ltransI, lthinI, fuseReductions);
1808
1809     if (fuseReductions) {
1810       BRepAlgoAPI_Fuse fuseInci (aResult, aReductionInci);
1811       if (!fuseInci.IsDone())
1812         StdFail_NotDone::Raise("Cannot fuse Te with incident reduction");
1813       aResult = fuseInci.Shape();
1814     }
1815     else {
1816       BRep_Builder B;
1817       TopoDS_Compound C;
1818       B.MakeCompound(C);
1819       B.Add(C, aResult);
1820       B.Add(C, aReductionInci);
1821       aResult = C;
1822     }
1823   }
1824
1825   // Get rid of extra compounds
1826   TopTools_ListOfShape listShapeRes;
1827   GEOMUtils::AddSimpleShapes(aResult, listShapeRes);
1828   aResult = listShapeRes.First(); // useful for the case "fuseReductions == true"
1829
1830   if (!fuseReductions && listShapeRes.Extent() > 1) {
1831     // Simplify T-Shape compound (get rid of sub-compounds) and glue duplicated faces
1832     BRep_Builder B;
1833     TopoDS_Compound C;
1834     B.MakeCompound(C);
1835
1836     TopTools_ListIteratorOfListOfShape itSub (listShapeRes);
1837     for (; itSub.More(); itSub.Next())
1838       B.Add(C, itSub.Value());
1839
1840     // GlueFaces
1841     aResult = GEOMImpl_GlueDriver::GlueFaces(C, Precision::Confusion(), Standard_True);
1842   }
1843
1844   return aResult;
1845 }
1846
1847 //=======================================================================
1848 //function : MakeThicknessReduction
1849 //purpose  : Static method. Create one thickness reduction element.
1850 //=======================================================================
1851 TopoDS_Shape GEOMImpl_IAdvancedOperations::MakeThicknessReduction (gp_Ax2 theAxes,
1852                                                                    const double R, const double W,
1853                                                                    const double Rthin, const double Wthin,
1854                                                                    const double Ltrans, const double Lthin,
1855                                                                    bool fuse)
1856 {
1857   double aTol = Precision::Confusion();
1858   if (Rthin < aTol || Wthin < aTol || Ltrans < aTol) {
1859     StdFail_NotDone::Raise("Cannot build thickness reduction: too small values");
1860   }
1861   bool isThinPart = (Lthin > aTol);
1862
1863   //     .
1864   //   W |\
1865   //     . \
1866   //   ^  \ '-----------------.
1867   //   |R  \|                 | Wthin
1868   //   |    '-----------------'
1869   //   v                        Rthin
1870   // --.--.--.--.--.--.--.--.--.--.--.--.--> theAxes.Direction()
1871   //     Ltrans     Lthin
1872
1873   double RExt = R + W;
1874   double RthinExt = Rthin + Wthin;
1875
1876   gp_Dir aNormal = theAxes.Direction();
1877   gp_Dir anXDir  = theAxes.XDirection();
1878   gp_Pnt aPntCyl (theAxes.Location().XYZ() + aNormal.XYZ()*Ltrans);
1879   gp_Ax2 anAxesCyl (aPntCyl, aNormal, anXDir);
1880
1881   // Build the transition part
1882   BRepPrimAPI_MakeCone ConeExt (theAxes, RExt, RthinExt, Ltrans);
1883   BRepPrimAPI_MakeCone ConeInt (theAxes, R, Rthin, Ltrans);
1884   ConeExt.Build();
1885   ConeInt.Build();
1886   if (!ConeExt.IsDone() || !ConeInt.IsDone())
1887     StdFail_NotDone::Raise("Cannot build cones of thickness reduction");
1888   BRepAlgoAPI_Cut cut1 (ConeExt.Shape(), ConeInt.Shape());
1889   if (!cut1.IsDone())
1890     StdFail_NotDone::Raise("Coudn't build transition part of thickness reduction");
1891   TopoDS_Shape aReduction = cut1.Shape();
1892
1893   // Build the thin part, if required
1894   TopoDS_Shape aThinPart;
1895   if (isThinPart) {
1896     BRepPrimAPI_MakeCylinder CExt (anAxesCyl, RthinExt, Lthin);
1897     BRepPrimAPI_MakeCylinder CInt (anAxesCyl, Rthin, Lthin);
1898     CExt.Build();
1899     CInt.Build();
1900     if (!CExt.IsDone() || !CInt.IsDone())
1901       StdFail_NotDone::Raise("Cannot build cylinders of thickness reduction");
1902     BRepAlgoAPI_Cut cut2 (CExt.Shape(), CInt.Shape());
1903     if (!cut2.IsDone())
1904       StdFail_NotDone::Raise("Coudn't build thin part of thickness reduction");
1905     aThinPart = cut2.Shape();
1906   }
1907
1908   // Join parts
1909   if (fuse) {
1910     if (isThinPart) {
1911       BRepAlgoAPI_Fuse fuse1 (aReduction, aThinPart);
1912       if (!fuse1.IsDone())
1913         StdFail_NotDone::Raise("Cannot fuse parts of thickness reduction");
1914       aReduction = fuse1.Shape();
1915     }
1916   }
1917   else {
1918     // Partition the reduction on blocks
1919     gp_Ax3 anAxesPln1 (aPntCyl, theAxes.XDirection(), aNormal);
1920     gp_Ax3 anAxesPln2 (aPntCyl, theAxes.YDirection(), aNormal);
1921     gp_Pln aPln1 (anAxesPln1);
1922     gp_Pln aPln2 (anAxesPln2);
1923     double aSize = Ltrans + Lthin + R + Rthin + Wthin; // to guarantee enough size in all directions
1924     TopoDS_Shape aTool1 = BRepBuilderAPI_MakeFace(aPln1, -aSize, +aSize, -aSize, +aSize).Shape();
1925     TopoDS_Shape aTool2 = BRepBuilderAPI_MakeFace(aPln2, -aSize, +aSize, -aSize, +aSize).Shape();
1926
1927     GEOMAlgo_Splitter PS;
1928     PS.AddArgument(aReduction);
1929     if (isThinPart)
1930       PS.AddArgument(aThinPart);
1931     PS.AddTool(aTool1);
1932     PS.AddTool(aTool2);
1933     PS.SetLimit(TopAbs_SOLID);
1934     PS.Perform();
1935
1936     aReduction = PS.Shape();
1937   }
1938
1939   return aReduction;
1940 }
1941
1942 //=============================================================================
1943 /*!
1944  *  MakePipeTShape
1945  *  \brief Create a T-shape object with specified caracteristics for the main and
1946  *         the incident pipes (radius, width, half-length).
1947  *         Center of the shape is (0,0,0). The main plane of the T-shape is XOY.
1948  *  \param theR1 Internal radius of main pipe
1949  *  \param theW1 Width of main pipe
1950  *  \param theL1 Half-length of main pipe
1951  *  \param theR2 Internal radius of incident pipe (R2 < R1)
1952  *  \param theW2 Width of incident pipe (R2+W2 < R1+W1)
1953  *  \param theL2 Half-length of incident pipe
1954  *  \param theHexMesh Boolean indicating if shape is prepared for hex mesh
1955  *  \return List of GEOM_Objects, containing the created shape and propagation groups.
1956  */
1957 //=============================================================================
1958 Handle(TColStd_HSequenceOfTransient)
1959   GEOMImpl_IAdvancedOperations::MakePipeTShape(double theR1, double theW1, double theL1,
1960                                                double theR2, double theW2, double theL2,
1961                                                double theRL, double theWL, double theLtransL, double theLthinL,
1962                                                double theRR, double theWR, double theLtransR, double theLthinR,
1963                                                double theRI, double theWI, double theLtransI, double theLthinI,
1964                                                bool theHexMesh)
1965 {
1966   MESSAGE("GEOMImpl_IAdvancedOperations::MakePipeTShape");
1967   SetErrorCode(KO);
1968   //Add a new object
1969   Handle(GEOM_Object) aShape = GetEngine()->AddObject(GetDocID(), GEOM_TSHAPE);
1970
1971   //Add a new shape function with parameters
1972   Handle(GEOM_Function) aFunction = aShape->AddFunction(GEOMImpl_PipeTShapeDriver::GetID(), TSHAPE_BASIC);
1973   if (aFunction.IsNull()) return NULL;
1974
1975   //Check if the function is set correctly
1976   if (aFunction->GetDriverGUID() != GEOMImpl_PipeTShapeDriver::GetID()) return NULL;
1977
1978   GEOMImpl_IPipeTShape aData (aFunction);
1979
1980   aData.SetR1(theR1);
1981   aData.SetW1(theW1);
1982   aData.SetL1(theL1);
1983   aData.SetR2(theR2);
1984   aData.SetW2(theW2);
1985   aData.SetL2(theL2);
1986   aData.SetHexMesh(theHexMesh);
1987
1988   bool isTRL = (theRL + theWL + theLtransL + theLthinL) > Precision::Confusion();
1989   bool isTRR = (theRR + theWR + theLtransR + theLthinR) > Precision::Confusion();
1990   bool isTRI = (theRI + theWI + theLtransI + theLthinI) > Precision::Confusion();
1991
1992   //Compute the resulting value
1993   try {
1994 #if OCC_VERSION_LARGE > 0x06010000
1995     OCC_CATCH_SIGNALS;
1996 #endif
1997     if (!GetSolver()->ComputeFunction(aFunction)) {
1998       SetErrorCode("TShape driver failed");
1999       return NULL;
2000     }
2001
2002     if (theHexMesh) {
2003       if (!MakePipeTShapePartition(aShape, theR1, theW1, theL1, theR2, theW2, theL2))
2004         return NULL;
2005       if (!MakePipeTShapeMirrorAndGlue(aShape, theR1, theW1, theL1, theR2, theW2, theL2))
2006         return NULL;
2007     }
2008
2009     if (isTRL || isTRR || isTRI) {
2010       // Add thickness reduction elements
2011       // at the three extremities: Left, Right and Incident
2012       TopoDS_Shape aResShape =
2013         MakePipeTShapeThicknessReduction(aShape->GetValue(), theR1, theW1, theL1, theR2, theW2, theL2,
2014                                          theRL, theWL, theLtransL, theLthinL,
2015                                          theRR, theWR, theLtransR, theLthinR,
2016                                          theRI, theWI, theLtransI, theLthinI,
2017                                          !theHexMesh);
2018       aFunction->SetValue(aResShape);
2019     }
2020   }
2021   catch (Standard_Failure) {
2022     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2023     SetErrorCode(aFail->GetMessageString());
2024     return NULL;
2025   }
2026
2027   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
2028   aSeq->Append(aShape);
2029
2030   try {
2031     if (theHexMesh) {
2032       // Get the groups
2033       if (!MakeGroups(aShape, TSHAPE_BASIC, theR1, theW1, theL1, theR2, theW2, theL2,
2034                       0., 0., 0., aSeq, gp_Trsf()))
2035         return NULL;
2036     }
2037
2038     // Get internal group.
2039     if (!MakeInternalGroup(aShape, theR1, theL1, theR2, theL2, theRL, theLtransL,
2040                            theRR, theLtransR, theRI, theLtransI,
2041                            aSeq, gp_Trsf())) {
2042       return NULL;
2043     }
2044   }
2045   catch (Standard_Failure) {
2046     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2047     SetErrorCode(aFail->GetMessageString());
2048     return NULL;
2049   }
2050
2051   //Make a Python command
2052   TCollection_AsciiString anEntry, aListRes("[");
2053   // Iterate over the sequence aSeq
2054   Standard_Integer aNbGroups = aSeq->Length();
2055   Standard_Integer i = 1;
2056   for (; i <= aNbGroups; i++) {
2057     Handle(Standard_Transient) anItem = aSeq->Value(i);
2058     if (anItem.IsNull()) continue;
2059     Handle(GEOM_Object) aGroup = Handle(GEOM_Object)::DownCast(anItem);
2060     if (aGroup.IsNull()) continue;
2061     //Make a Python command
2062     TDF_Tool::Entry(aGroup->GetEntry(), anEntry);
2063     aListRes += anEntry + ", ";
2064   }
2065   aListRes.Trunc(aListRes.Length() - 2);
2066
2067   GEOM::TPythonDump pd (aFunction);
2068
2069   pd << aListRes.ToCString() << "] = geompy.MakePipeTShape("
2070      << theR1 << ", " << theW1 << ", " << theL1 << ", "
2071      << theR2 << ", " << theW2 << ", " << theL2 << ", "
2072      << theHexMesh;
2073
2074   // thickness reduction
2075   if (isTRL)
2076     pd << ", theRL=" << theRL << ", theWL=" << theWL
2077        << ", theLtransL=" << theLtransL << ", theLthinL=" << theLthinL;
2078   if (isTRR)
2079     pd << ", theRR=" << theRR << ", theWR=" << theWR
2080        << ", theLtransR=" << theLtransR << ", theLthinR=" << theLthinR;
2081   if (isTRI)
2082     pd << ", theRI=" << theRI << ", theWI=" << theWI
2083        << ", theLtransI=" << theLtransI << ", theLthinI=" << theLthinI;
2084
2085   pd << ")";
2086
2087   SetErrorCode(OK);
2088
2089   return aSeq;
2090 }
2091
2092 //=============================================================================
2093 /*!
2094  *  MakePipeTShapeWithPosition
2095  *  Create a T-shape object with specified caracteristics for the main and
2096  *  the incident pipes (radius, width, half-length).
2097  *  The extremities of the main pipe are located on junctions points P1 and P2.
2098  *  The extremity of the incident pipe is located on junction point P3.
2099  *  \param theR1 Internal radius of main pipe
2100  *  \param theW1 Width of main pipe
2101  *  \param theL1 Half-length of main pipe
2102  *  \param theR2 Internal radius of incident pipe (R2 < R1)
2103  *  \param theW2 Width of incident pipe (R2+W2 < R1+W1)
2104  *  \param theL2 Half-length of incident pipe
2105  *  \param theHexMesh Boolean indicating if shape is prepared for hex mesh
2106  *  \param theP1 1st junction point of main pipe
2107  *  \param theP2 2nd junction point of main pipe
2108  *  \param theP3 Junction point of incident pipe
2109  *  \return List of GEOM_Objects, containing the created shape and propagation groups..
2110  */
2111 //=============================================================================
2112 Handle(TColStd_HSequenceOfTransient)
2113 GEOMImpl_IAdvancedOperations::MakePipeTShapeWithPosition
2114                              (double theR1, double theW1, double theL1,
2115                               double theR2, double theW2, double theL2,
2116                               double theRL, double theWL, double theLtransL, double theLthinL,
2117                               double theRR, double theWR, double theLtransR, double theLthinR,
2118                               double theRI, double theWI, double theLtransI, double theLthinI,
2119                               bool theHexMesh,
2120                               Handle(GEOM_Object) theP1,
2121                               Handle(GEOM_Object) theP2,
2122                               Handle(GEOM_Object) theP3)
2123 {
2124   SetErrorCode(KO);
2125   //Add a new object
2126   Handle(GEOM_Object) aShape = GetEngine()->AddObject(GetDocID(), GEOM_TSHAPE);
2127   /////////////////
2128   // TSHAPE CODE
2129   /////////////////
2130   //Add a new shape function with parameters
2131   Handle(GEOM_Function) aFunction = aShape->AddFunction(GEOMImpl_PipeTShapeDriver::GetID(), TSHAPE_BASIC);
2132   if (aFunction.IsNull()) return NULL;
2133
2134   //Check if the function is set correctly
2135   if (aFunction->GetDriverGUID() != GEOMImpl_PipeTShapeDriver::GetID()) return NULL;
2136
2137   // Check new position
2138   if (!CheckCompatiblePosition(theL1, theL2, theP1, theP2, theP3, 0.01)) {
2139     return NULL;
2140   }
2141
2142   GEOMImpl_IPipeTShape aData(aFunction);
2143
2144   aData.SetR1(theR1);
2145   aData.SetW1(theW1);
2146   aData.SetL1(theL1);
2147   aData.SetR2(theR2);
2148   aData.SetW2(theW2);
2149   aData.SetL2(theL2);
2150   aData.SetHexMesh(theHexMesh);
2151
2152   bool isTRL = (theRL + theWL + theLtransL + theLthinL) > Precision::Confusion();
2153   bool isTRR = (theRR + theWR + theLtransR + theLthinR) > Precision::Confusion();
2154   bool isTRI = (theRI + theWI + theLtransI + theLthinI) > Precision::Confusion();
2155
2156   //Compute the resulting value
2157   try {
2158 #if OCC_VERSION_LARGE > 0x06010000
2159     OCC_CATCH_SIGNALS;
2160 #endif
2161     if (!GetSolver()->ComputeFunction(aFunction)) {
2162       SetErrorCode("TShape driver failed");
2163       return NULL;
2164     }
2165
2166     if (theHexMesh) {
2167       if (!MakePipeTShapePartition(aShape, theR1, theW1, theL1, theR2, theW2, theL2))
2168         return NULL;
2169       if (!MakePipeTShapeMirrorAndGlue(aShape, theR1, theW1, theL1, theR2, theW2, theL2))
2170         return NULL;
2171     }
2172
2173     if (isTRL || isTRR || isTRI) {
2174       // Add thickness reduction elements
2175       // at the three extremities: Left, Right and Incident
2176       TopoDS_Shape aResShape =
2177         MakePipeTShapeThicknessReduction(aShape->GetValue(), theR1, theW1, theL1, theR2, theW2, theL2,
2178                                          theRL, theWL, theLtransL, theLthinL,
2179                                          theRR, theWR, theLtransR, theLthinR,
2180                                          theRI, theWI, theLtransI, theLthinI,
2181                                          !theHexMesh);
2182       aFunction->SetValue(aResShape);
2183     }
2184   }
2185   catch (Standard_Failure) {
2186     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2187     SetErrorCode(aFail->GetMessageString());
2188     return NULL;
2189   }
2190
2191   TopoDS_Shape Te = aShape->GetValue();
2192
2193   // Set Position
2194   gp_Trsf aTrsf = GetPositionTrsf(theL1, theL2, theP1, theP2, theP3);
2195   BRepBuilderAPI_Transform aTransformation(Te, aTrsf, Standard_False);
2196   TopoDS_Shape aTrsf_Shape = aTransformation.Shape();
2197   aFunction->SetValue(aTrsf_Shape);
2198
2199   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
2200   aSeq->Append(aShape);
2201
2202   try {
2203     if (theHexMesh) {
2204       // Get the groups
2205       if (!MakeGroups(aShape,TSHAPE_BASIC, theR1, theW1, theL1, theR2, theW2, theL2,
2206                       0., 0., 0., aSeq, aTrsf)) {
2207         return NULL;
2208       }
2209     }
2210
2211     // Get internal group.
2212     if (!MakeInternalGroup(aShape, theR1, theL1, theR2, theL2, theRL, theLtransL,
2213                            theRR, theLtransR, theRI, theLtransI,
2214                            aSeq, aTrsf)) {
2215       return NULL;
2216     }
2217   }
2218   catch (Standard_Failure) {
2219     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2220     SetErrorCode(aFail->GetMessageString());
2221     return NULL;
2222   }
2223
2224   //Make a Python command
2225   TCollection_AsciiString anEntry, aListRes("[");
2226   // Iterate over the sequence aSeq
2227   Standard_Integer aNbGroups = aSeq->Length();
2228   Standard_Integer i = 1;
2229   for (; i <= aNbGroups; i++) {
2230     Handle(Standard_Transient) anItem = aSeq->Value(i);
2231     if (anItem.IsNull()) continue;
2232     Handle(GEOM_Object) aGroup = Handle(GEOM_Object)::DownCast(anItem);
2233     if (aGroup.IsNull()) continue;
2234     //Make a Python command
2235     TDF_Tool::Entry(aGroup->GetEntry(), anEntry);
2236     aListRes += anEntry + ", ";
2237   }
2238   aListRes.Trunc(aListRes.Length() - 2);
2239
2240   GEOM::TPythonDump pd (aFunction);
2241
2242   pd << aListRes.ToCString() << "] = geompy.MakePipeTShape("
2243      << theR1 << ", " << theW1 << ", " << theL1 << ", "
2244      << theR2 << ", " << theW2 << ", " << theL2 << ", "
2245      << theHexMesh << ", " << theP1 << ", " << theP2 << ", " << theP3;
2246
2247   // thickness reduction
2248   if (isTRL)
2249     pd << ", theRL=" << theRL << ", theWL=" << theWL
2250        << ", theLtransL=" << theLtransL << ", theLthinL=" << theLthinL;
2251   if (isTRR)
2252     pd << ", theRR=" << theRR << ", theWR=" << theWR
2253        << ", theLtransR=" << theLtransR << ", theLthinR=" << theLthinR;
2254   if (isTRI)
2255     pd << ", theRI=" << theRI << ", theWI=" << theWI
2256        << ", theLtransI=" << theLtransI << ", theLthinI=" << theLthinI;
2257
2258   pd << ")";
2259
2260   SetErrorCode(OK);
2261
2262   return aSeq;
2263 }
2264
2265 //=============================================================================
2266 /*!
2267  *  MakePipeTShapeChamfer
2268  *  Create a T-shape object with specified caracteristics for the main and
2269  *  the incident pipes (radius, width, half-length). A chamfer is created
2270  *  on the junction of the pipes.
2271  *  Center of the shape is (0,0,0). The main plane of the T-shape is XOY.
2272  *  \param theR1 Internal radius of main pipe
2273  *  \param theW1 Width of main pipe
2274  *  \param theL1 Half-length of main pipe
2275  *  \param theR2 Internal radius of incident pipe (R2 < R1)
2276  *  \param theW2 Width of incident pipe (R2+W2 < R1+W1)
2277  *  \param theL2 Half-length of incident pipe
2278  *  \param theH Height of chamfer.
2279  *  \param theW Width of chamfer.
2280  *  \param theHexMesh Boolean indicating if shape is prepared for hex mesh
2281  *  \return List of GEOM_Objects, containing the created shape and propagation groups.
2282  */
2283 //=============================================================================
2284 Handle(TColStd_HSequenceOfTransient)
2285 GEOMImpl_IAdvancedOperations::MakePipeTShapeChamfer
2286                              (double theR1, double theW1, double theL1,
2287                               double theR2, double theW2, double theL2,
2288                               double theRL, double theWL, double theLtransL, double theLthinL,
2289                               double theRR, double theWR, double theLtransR, double theLthinR,
2290                               double theRI, double theWI, double theLtransI, double theLthinI,
2291                               double theH, double theW,
2292                               bool theHexMesh)
2293 {
2294   SetErrorCode(KO);
2295   //Add a new object
2296   Handle(GEOM_Object) aShape = GetEngine()->AddObject(GetDocID(), GEOM_TSHAPE);
2297   //Add a new shape function with parameters
2298   Handle(GEOM_Function) aFunction = aShape->AddFunction(GEOMImpl_PipeTShapeDriver::GetID(), TSHAPE_CHAMFER);
2299   if (aFunction.IsNull()) return NULL;
2300
2301   //Check if the function is set correctly
2302   if (aFunction->GetDriverGUID() != GEOMImpl_PipeTShapeDriver::GetID()) return NULL;
2303
2304   GEOMImpl_IPipeTShape aData(aFunction);
2305
2306   aData.SetR1(theR1);
2307   aData.SetW1(theW1);
2308   aData.SetL1(theL1);
2309   aData.SetR2(theR2);
2310   aData.SetW2(theW2);
2311   aData.SetL2(theL2);
2312   aData.SetH(theH);
2313   aData.SetW(theW);
2314   aData.SetHexMesh(theHexMesh);
2315
2316   bool isTRL = (theRL + theWL + theLtransL + theLthinL) > Precision::Confusion();
2317   bool isTRR = (theRR + theWR + theLtransR + theLthinR) > Precision::Confusion();
2318   bool isTRI = (theRI + theWI + theLtransI + theLthinI) > Precision::Confusion();
2319
2320   //Compute the resulting value
2321   try {
2322 #if OCC_VERSION_LARGE > 0x06010000
2323     OCC_CATCH_SIGNALS;
2324 #endif
2325     if (!GetSolver()->ComputeFunction(aFunction)) {
2326       SetErrorCode("TShape driver failed");
2327       return NULL;
2328     }
2329   }
2330   catch (Standard_Failure) {
2331     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2332     SetErrorCode(aFail->GetMessageString());
2333     return NULL;
2334   }
2335
2336   // BEGIN of chamfer
2337   TopoDS_Shape aShapeShape = aShape->GetValue();
2338   TopTools_IndexedMapOfShape anEdgesIndices;
2339   TopExp::MapShapes(aShapeShape, anEdgesIndices);
2340   // Common edges on external cylinders
2341   Handle(GEOM_Object) box_e;
2342   if (theHexMesh) {
2343     box_e = my3DPrimOperations->MakeBoxDXDYDZ(theR2+theW2, theR2+theW2, theR1+theW1);
2344   }
2345   else {
2346     box_e = my3DPrimOperations->MakeBoxDXDYDZ(2*(theR2+theW2), 2*(theR2+theW2), theR1+theW1);
2347   }
2348   box_e->GetLastFunction()->SetDescription("");
2349   box_e = myTransformOperations->TranslateDXDYDZ(box_e, -(theR2+theW2), -(theR2+theW2), 0);
2350   box_e->GetLastFunction()->SetDescription("");
2351
2352   Handle(TColStd_HSequenceOfInteger) edges_e =
2353     myShapesOperations->GetShapesOnBoxIDs(box_e, aShape, TopAbs_EDGE, GEOMAlgo_ST_IN);
2354   box_e->GetLastFunction()->SetDescription("");
2355
2356   if (edges_e.IsNull() || edges_e->Length() == 0) {
2357     SetErrorCode("External edges not found");
2358     return NULL;
2359   }
2360   int nbEdgesInChamfer = 0;
2361   std::list<int> theEdges;
2362   for (int i=1; i<=edges_e->Length();i++) {
2363     int edgeID = edges_e->Value(i);
2364     TopoDS_Shape theEdge = anEdgesIndices.FindKey(edgeID);
2365     TopExp_Explorer Ex(theEdge,TopAbs_VERTEX);
2366     int iv=0;
2367     while (Ex.More()) {
2368       iv ++;
2369       gp_Pnt aPt = BRep_Tool::Pnt(TopoDS::Vertex(Ex.Current()));
2370       if (Abs(aPt.Z() - (theR1+theW1)) <= Precision::Confusion()) {
2371         nbEdgesInChamfer ++;
2372         theEdges.push_back(edgeID);
2373       }
2374       Ex.Next();
2375     }
2376     if (theHexMesh && nbEdgesInChamfer == 1)
2377       break;
2378   }
2379   Handle(GEOM_Object) aChamfer;
2380   try {
2381     aChamfer = myLocalOperations->MakeChamferEdges(aShape, theW, theH, theEdges);
2382   }
2383   catch (Standard_Failure) {
2384     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2385     SetErrorCode(aFail->GetMessageString());
2386     return NULL;
2387   }
2388   if (aChamfer.IsNull()) {
2389     SetErrorCode("Chamfer can not be computed on the given shape with the given parameters");
2390     return NULL;
2391   }
2392   aChamfer->GetLastFunction()->SetDescription("");
2393
2394   TopoDS_Shape aChamferShape = aChamfer->GetValue();
2395   aFunction->SetValue(aChamferShape);
2396   // END of chamfer
2397
2398   if (theHexMesh) {
2399     if (!MakePipeTShapePartition(aShape, theR1, theW1, theL1, theR2, theW2, theL2, theH, theW, 0, false))
2400       return NULL;
2401     if (!MakePipeTShapeMirrorAndGlue(aShape, theR1, theW1, theL1, theR2, theW2, theL2))
2402       return NULL;
2403   }
2404
2405   // Add thickness reduction elements
2406   // at the three extremities: Left, Right and Incident
2407   try {
2408 #if OCC_VERSION_LARGE > 0x06010000
2409     OCC_CATCH_SIGNALS;
2410 #endif
2411     if (isTRL || isTRR || isTRI) {
2412       TopoDS_Shape aResShape =
2413         MakePipeTShapeThicknessReduction(aShape->GetValue(), theR1, theW1, theL1, theR2, theW2, theL2,
2414                                          theRL, theWL, theLtransL, theLthinL,
2415                                          theRR, theWR, theLtransR, theLthinR,
2416                                          theRI, theWI, theLtransI, theLthinI,
2417                                          !theHexMesh);
2418       aFunction->SetValue(aResShape);
2419     }
2420   }
2421   catch (Standard_Failure) {
2422     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2423     SetErrorCode(aFail->GetMessageString());
2424     return NULL;
2425   }
2426
2427   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
2428   aSeq->Append(aShape);
2429
2430   try {
2431     if (theHexMesh) {
2432       // Get the groups
2433       if (!MakeGroups(aShape, TSHAPE_CHAMFER, theR1, theW1, theL1, theR2, theW2, theL2,
2434                       theH, theW, 0., aSeq, gp_Trsf()))
2435         return NULL;
2436     }
2437
2438     // Get internal group.
2439     if (!MakeInternalGroup(aShape, theR1, theL1, theR2, theL2, theRL, theLtransL,
2440                            theRR, theLtransR, theRI, theLtransI,
2441                            aSeq, gp_Trsf())) {
2442       return NULL;
2443     }
2444   }
2445   catch (Standard_Failure) {
2446     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2447     SetErrorCode(aFail->GetMessageString());
2448     return NULL;
2449   }
2450
2451   //Make a Python command
2452   TCollection_AsciiString anEntry, aListRes("[");
2453   // Iterate over the sequence aSeq
2454   Standard_Integer aNbGroups = aSeq->Length();
2455   Standard_Integer i = 1;
2456   for (; i <= aNbGroups; i++) {
2457     Handle(Standard_Transient) anItem = aSeq->Value(i);
2458     if (anItem.IsNull()) continue;
2459     Handle(GEOM_Object) aGroup = Handle(GEOM_Object)::DownCast(anItem);
2460     if (aGroup.IsNull()) continue;
2461     //Make a Python command
2462     TDF_Tool::Entry(aGroup->GetEntry(), anEntry);
2463     aListRes += anEntry + ", ";
2464   }
2465   aListRes.Trunc(aListRes.Length() - 2);
2466
2467   GEOM::TPythonDump pd (aFunction);
2468
2469   pd << aListRes.ToCString() << "] = geompy.MakePipeTShapeChamfer("
2470      << theR1 << ", " << theW1 << ", " << theL1 << ", "
2471      << theR2 << ", " << theW2 << ", " << theL2 << ", "
2472      << theH << ", " << theW << ", " << theHexMesh;
2473
2474   // thickness reduction
2475   if (isTRL)
2476     pd << ", theRL=" << theRL << ", theWL=" << theWL
2477        << ", theLtransL=" << theLtransL << ", theLthinL=" << theLthinL;
2478   if (isTRR)
2479     pd << ", theRR=" << theRR << ", theWR=" << theWR
2480        << ", theLtransR=" << theLtransR << ", theLthinR=" << theLthinR;
2481   if (isTRI)
2482     pd << ", theRI=" << theRI << ", theWI=" << theWI
2483        << ", theLtransI=" << theLtransI << ", theLthinI=" << theLthinI;
2484
2485   pd << ")";
2486
2487   SetErrorCode(OK);
2488
2489   return aSeq;
2490 }
2491
2492 //=============================================================================
2493 /*!
2494  *  MakePipeTShapeChamferWithPosition
2495  *  Create a T-shape object with specified caracteristics for the main and
2496  *  the incident pipes (radius, width, half-length). A chamfer is created
2497  *  on the junction of the pipes.
2498  *  The extremities of the main pipe are located on junctions points P1 and P2.
2499  *  The extremity of the incident pipe is located on junction point P3.
2500  *  \param theR1 Internal radius of main pipe
2501  *  \param theW1 Width of main pipe
2502  *  \param theL1 Half-length of main pipe
2503  *  \param theR2 Internal radius of incident pipe (R2 < R1)
2504  *  \param theW2 Width of incident pipe (R2+W2 < R1+W1)
2505  *  \param theL2 Half-length of incident pipe
2506  *  \param theH Height of chamfer.
2507  *  \param theW Width of chamfer.
2508  *  \param theHexMesh Boolean indicating if shape is prepared for hex mesh
2509  *  \param theP1 1st junction point of main pipe
2510  *  \param theP2 2nd junction point of main pipe
2511  *  \param theP3 Junction point of incident pipe
2512  *  \return List of GEOM_Objects, containing the created shape and propagation groups.
2513  */
2514 //=============================================================================
2515 Handle(TColStd_HSequenceOfTransient)
2516 GEOMImpl_IAdvancedOperations::MakePipeTShapeChamferWithPosition
2517                              (double theR1, double theW1, double theL1,
2518                               double theR2, double theW2, double theL2,
2519                               double theRL, double theWL, double theLtransL, double theLthinL,
2520                               double theRR, double theWR, double theLtransR, double theLthinR,
2521                               double theRI, double theWI, double theLtransI, double theLthinI,
2522                               double theH, double theW,
2523                               bool theHexMesh,
2524                               Handle(GEOM_Object) theP1,
2525                               Handle(GEOM_Object) theP2,
2526                               Handle(GEOM_Object) theP3)
2527 {
2528   SetErrorCode(KO);
2529   //Add a new object
2530   Handle(GEOM_Object) aShape = GetEngine()->AddObject(GetDocID(), GEOM_TSHAPE);
2531   //Add a new shape function with parameters
2532   Handle(GEOM_Function) aFunction = aShape->AddFunction(GEOMImpl_PipeTShapeDriver::GetID(), TSHAPE_CHAMFER);
2533   if (aFunction.IsNull()) return NULL;
2534
2535   //Check if the function is set correctly
2536   if (aFunction->GetDriverGUID() != GEOMImpl_PipeTShapeDriver::GetID()) return NULL;
2537
2538   // Check new position
2539   if (!CheckCompatiblePosition(theL1, theL2, theP1, theP2, theP3, 0.01)) {
2540     return NULL;
2541   }
2542
2543   GEOMImpl_IPipeTShape aData(aFunction);
2544
2545   aData.SetR1(theR1);
2546   aData.SetW1(theW1);
2547   aData.SetL1(theL1);
2548   aData.SetR2(theR2);
2549   aData.SetW2(theW2);
2550   aData.SetL2(theL2);
2551   aData.SetH(theH);
2552   aData.SetW(theW);
2553   aData.SetHexMesh(theHexMesh);
2554
2555   bool isTRL = (theRL + theWL + theLtransL + theLthinL) > Precision::Confusion();
2556   bool isTRR = (theRR + theWR + theLtransR + theLthinR) > Precision::Confusion();
2557   bool isTRI = (theRI + theWI + theLtransI + theLthinI) > Precision::Confusion();
2558
2559   //Compute the resulting value
2560   try {
2561 #if OCC_VERSION_LARGE > 0x06010000
2562     OCC_CATCH_SIGNALS;
2563 #endif
2564     if (!GetSolver()->ComputeFunction(aFunction)) {
2565       SetErrorCode("TShape driver failed");
2566       return NULL;
2567     }
2568   }
2569   catch (Standard_Failure) {
2570     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2571     SetErrorCode(aFail->GetMessageString());
2572     return NULL;
2573   }
2574
2575   // BEGIN of chamfer
2576   TopoDS_Shape aShapeShape = aShape->GetValue();
2577   TopTools_IndexedMapOfShape anEdgesIndices;
2578   TopExp::MapShapes(aShapeShape, anEdgesIndices);
2579   // Common edges on external cylinders
2580   Handle(GEOM_Object) box_e;
2581   if (theHexMesh) {
2582     box_e = my3DPrimOperations->MakeBoxDXDYDZ(theR2+theW2, theR2+theW2, theR1+theW1);
2583   }
2584   else {
2585     box_e = my3DPrimOperations->MakeBoxDXDYDZ(2*(theR2+theW2), 2*(theR2+theW2), theR1+theW1);
2586   }
2587   box_e->GetLastFunction()->SetDescription("");
2588   box_e = myTransformOperations->TranslateDXDYDZ(box_e, -(theR2+theW2), -(theR2+theW2), 0);
2589   box_e->GetLastFunction()->SetDescription("");
2590
2591   Handle(TColStd_HSequenceOfInteger) edges_e =
2592     myShapesOperations->GetShapesOnBoxIDs(box_e, aShape, TopAbs_EDGE, GEOMAlgo_ST_IN);
2593   box_e->GetLastFunction()->SetDescription("");
2594
2595   if (edges_e.IsNull() || edges_e->Length() == 0) {
2596     SetErrorCode("External edges not found");
2597     return NULL;
2598   }
2599   int nbEdgesInChamfer = 0;
2600   std::list<int> theEdges;
2601   for (int i=1; i<=edges_e->Length();i++) {
2602     int edgeID = edges_e->Value(i);
2603     TopoDS_Shape theEdge = anEdgesIndices.FindKey(edgeID);
2604     TopExp_Explorer Ex(theEdge,TopAbs_VERTEX);
2605     while (Ex.More()) {
2606       gp_Pnt aPt = BRep_Tool::Pnt(TopoDS::Vertex(Ex.Current()));
2607       if (Abs(aPt.Z() - (theR1+theW1)) <= Precision::Confusion()) {
2608         nbEdgesInChamfer ++;
2609         theEdges.push_back(edgeID);
2610       }
2611       Ex.Next();
2612     }
2613     if (theHexMesh && nbEdgesInChamfer == 1)
2614       break;
2615   }
2616   Handle(GEOM_Object) aChamfer;
2617   try {
2618     aChamfer = myLocalOperations->MakeChamferEdges(aShape, theW, theH, theEdges);
2619   }
2620   catch (Standard_Failure) {
2621     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2622     SetErrorCode(aFail->GetMessageString());
2623     return NULL;
2624   }
2625   if (aChamfer.IsNull()) {
2626     SetErrorCode("Chamfer can not be computed on the given shape with the given parameters");
2627     return NULL;
2628   }
2629   aChamfer->GetLastFunction()->SetDescription("");
2630
2631   TopoDS_Shape aChamferShape = aChamfer->GetValue();
2632   aFunction->SetValue(aChamferShape);
2633   // END of chamfer
2634
2635   if (theHexMesh) {
2636     if (!MakePipeTShapePartition(aShape, theR1, theW1, theL1, theR2, theW2, theL2, theH, theW, 0, false))
2637       return NULL;
2638     if (!MakePipeTShapeMirrorAndGlue(aShape, theR1, theW1, theL1, theR2, theW2, theL2))
2639       return NULL;
2640   }
2641
2642   // Add thickness reduction elements
2643   // at the three extremities: Left, Right and Incident
2644   try {
2645 #if OCC_VERSION_LARGE > 0x06010000
2646     OCC_CATCH_SIGNALS;
2647 #endif
2648     if (isTRL || isTRR || isTRI) {
2649       TopoDS_Shape aResShape =
2650         MakePipeTShapeThicknessReduction(aShape->GetValue(), theR1, theW1, theL1, theR2, theW2, theL2,
2651                                          theRL, theWL, theLtransL, theLthinL,
2652                                          theRR, theWR, theLtransR, theLthinR,
2653                                          theRI, theWI, theLtransI, theLthinI,
2654                                          !theHexMesh);
2655       aFunction->SetValue(aResShape);
2656     }
2657   }
2658   catch (Standard_Failure) {
2659     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2660     SetErrorCode(aFail->GetMessageString());
2661     return NULL;
2662   }
2663
2664   // Set Position
2665   gp_Trsf aTrsf = GetPositionTrsf(theL1, theL2, theP1, theP2, theP3);
2666   BRepBuilderAPI_Transform aTransformation (aShape->GetValue(), aTrsf, Standard_False);
2667   TopoDS_Shape aTrsf_Shape = aTransformation.Shape();
2668   aFunction->SetValue(aTrsf_Shape);
2669
2670   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
2671   aSeq->Append(aShape);
2672
2673   try {
2674     if (theHexMesh) {
2675       // Get the groups
2676       if (!MakeGroups(aShape, TSHAPE_CHAMFER, theR1, theW1, theL1, theR2, theW2, theL2,
2677                       theH, theW, 0., aSeq, aTrsf))
2678         return NULL;
2679     }
2680
2681     // Get internal group.
2682     if (!MakeInternalGroup(aShape, theR1, theL1, theR2, theL2, theRL, theLtransL,
2683                            theRR, theLtransR, theRI, theLtransI,
2684                            aSeq, aTrsf)) {
2685       return NULL;
2686     }
2687   }
2688   catch (Standard_Failure) {
2689     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2690     SetErrorCode(aFail->GetMessageString());
2691     return NULL;
2692   }
2693
2694   //Make a Python command
2695   TCollection_AsciiString anEntry, aListRes("[");
2696   // Iterate over the sequence aSeq
2697   Standard_Integer aNbGroups = aSeq->Length();
2698   Standard_Integer i = 1;
2699   for (; i <= aNbGroups; i++) {
2700     Handle(Standard_Transient) anItem = aSeq->Value(i);
2701     if (anItem.IsNull()) continue;
2702     Handle(GEOM_Object) aGroup = Handle(GEOM_Object)::DownCast(anItem);
2703     if (aGroup.IsNull()) continue;
2704     //Make a Python command
2705     TDF_Tool::Entry(aGroup->GetEntry(), anEntry);
2706     aListRes += anEntry + ", ";
2707   }
2708   aListRes.Trunc(aListRes.Length() - 2);
2709
2710   GEOM::TPythonDump pd (aFunction);
2711
2712   pd << aListRes.ToCString() << "] = geompy.MakePipeTShapeChamfer("
2713      << theR1 << ", " << theW1 << ", " << theL1 << ", "
2714      << theR2 << ", " << theW2 << ", " << theL2 << ", "
2715      << theH << ", " << theW << ", " << theHexMesh << ", "
2716      << theP1 << ", " << theP2 << ", " << theP3;
2717
2718   // thickness reduction
2719   if (isTRL)
2720     pd << ", theRL=" << theRL << ", theWL=" << theWL
2721        << ", theLtransL=" << theLtransL << ", theLthinL=" << theLthinL;
2722   if (isTRR)
2723     pd << ", theRR=" << theRR << ", theWR=" << theWR
2724        << ", theLtransR=" << theLtransR << ", theLthinR=" << theLthinR;
2725   if (isTRI)
2726     pd << ", theRI=" << theRI << ", theWI=" << theWI
2727        << ", theLtransI=" << theLtransI << ", theLthinI=" << theLthinI;
2728
2729   pd << ")";
2730
2731   SetErrorCode(OK);
2732
2733   return aSeq;
2734 }
2735
2736 //=============================================================================
2737 /*!
2738  *  MakePipeTShapeFillet
2739  *  Create a T-shape object with specified caracteristics for the main and
2740  *  the incident pipes (radius, width, half-length). A fillet is created
2741  *  on the junction of the pipes.
2742  *  Center of the shape is (0,0,0). The main plane of the T-shape is XOY.
2743  *  \param theR1 Internal radius of main pipe
2744  *  \param theW1 Width of main pipe
2745  *  \param theL1 Half-length of main pipe
2746  *  \param theR2 Internal radius of incident pipe (R2 < R1)
2747  *  \param theW2 Width of incident pipe (R2+W2 < R1+W1)
2748  *  \param theL2 Half-length of incident pipe
2749  *  \param theRF Radius of curvature of fillet.
2750  *  \param theHexMesh Boolean indicating if shape is prepared for hex mesh
2751  *  \return List of GEOM_Objects, containing the created shape and propagation groups.
2752  */
2753 //=============================================================================
2754 Handle(TColStd_HSequenceOfTransient)
2755 GEOMImpl_IAdvancedOperations::MakePipeTShapeFillet
2756                              (double theR1, double theW1, double theL1,
2757                               double theR2, double theW2, double theL2,
2758                               double theRL, double theWL, double theLtransL, double theLthinL,
2759                               double theRR, double theWR, double theLtransR, double theLthinR,
2760                               double theRI, double theWI, double theLtransI, double theLthinI,
2761                               double theRF, bool theHexMesh)
2762 {
2763   SetErrorCode(KO);
2764   //Add a new object
2765   Handle(GEOM_Object) aShape = GetEngine()->AddObject(GetDocID(), GEOM_TSHAPE);
2766   //Add a new shape function with parameters
2767   Handle(GEOM_Function) aFunction = aShape->AddFunction(GEOMImpl_PipeTShapeDriver::GetID(), TSHAPE_FILLET);
2768   if (aFunction.IsNull()) return NULL;
2769
2770   //Check if the function is set correctly
2771   if (aFunction->GetDriverGUID() != GEOMImpl_PipeTShapeDriver::GetID()) return NULL;
2772
2773   GEOMImpl_IPipeTShape aData(aFunction);
2774
2775   aData.SetR1(theR1);
2776   aData.SetW1(theW1);
2777   aData.SetL1(theL1);
2778   aData.SetR2(theR2);
2779   aData.SetW2(theW2);
2780   aData.SetL2(theL2);
2781   aData.SetRF(theRF);
2782   aData.SetHexMesh(theHexMesh);
2783
2784   bool isTRL = (theRL + theWL + theLtransL + theLthinL) > Precision::Confusion();
2785   bool isTRR = (theRR + theWR + theLtransR + theLthinR) > Precision::Confusion();
2786   bool isTRI = (theRI + theWI + theLtransI + theLthinI) > Precision::Confusion();
2787
2788   //Compute the resulting value
2789   try {
2790 #if OCC_VERSION_LARGE > 0x06010000
2791     OCC_CATCH_SIGNALS;
2792 #endif
2793     if (!GetSolver()->ComputeFunction(aFunction)) {
2794       SetErrorCode("TShape driver failed");
2795       return NULL;
2796     }
2797   }
2798   catch (Standard_Failure) {
2799     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2800     SetErrorCode(aFail->GetMessageString());
2801     return NULL;
2802   }
2803
2804   // BEGIN of fillet
2805   TopoDS_Shape aShapeShape = aShape->GetValue();
2806   TopTools_IndexedMapOfShape anEdgesIndices;
2807   TopExp::MapShapes(aShapeShape, anEdgesIndices);
2808   // Common edges on external cylinders
2809   Handle(GEOM_Object) box_e;
2810   if (theHexMesh) {
2811     box_e = my3DPrimOperations->MakeBoxDXDYDZ(theR2+theW2, theR2+theW2, theR1+theW1);
2812   }
2813   else {
2814     box_e = my3DPrimOperations->MakeBoxDXDYDZ(2*(theR2+theW2), 2*(theR2+theW2), theR1+theW1);
2815   }
2816   box_e->GetLastFunction()->SetDescription("");
2817   box_e = myTransformOperations->TranslateDXDYDZ(box_e, -(theR2+theW2), -(theR2+theW2), 0);
2818   box_e->GetLastFunction()->SetDescription("");
2819
2820   Handle(TColStd_HSequenceOfInteger) edges_e =
2821     myShapesOperations->GetShapesOnBoxIDs(box_e, aShape, TopAbs_EDGE, GEOMAlgo_ST_IN);
2822   box_e->GetLastFunction()->SetDescription("");
2823
2824   if (edges_e.IsNull() || edges_e->Length() == 0) {
2825     SetErrorCode("External edges not found");
2826     return NULL;
2827   }
2828   int nbEdgesInFillet = 0;
2829   std::list<int> theEdges;
2830   for (int i=1; i<=edges_e->Length();i++) {
2831     int edgeID = edges_e->Value(i);
2832     TopoDS_Shape theEdge = anEdgesIndices.FindKey(edgeID);
2833     TopExp_Explorer Ex(theEdge,TopAbs_VERTEX);
2834     while (Ex.More()) {
2835       gp_Pnt aPt = BRep_Tool::Pnt(TopoDS::Vertex(Ex.Current()));
2836       if (Abs(aPt.Z() - (theR1+theW1)) <= Precision::Confusion()) {
2837         nbEdgesInFillet ++;
2838         theEdges.push_back(edgeID);
2839       }
2840       Ex.Next();
2841     }
2842     if (theHexMesh && nbEdgesInFillet == 1)
2843       break;
2844   }
2845
2846   Handle(GEOM_Object) aFillet;
2847   try {
2848     aFillet = myLocalOperations->MakeFilletEdges(aShape, theRF, theEdges);
2849   }
2850   catch (Standard_Failure) {
2851     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2852     SetErrorCode(aFail->GetMessageString());
2853     return NULL;
2854   }
2855   if (aFillet.IsNull()) {
2856     //SetErrorCode("Fillet can not be computed on the given shape with the given parameters");
2857     SetErrorCode(myLocalOperations->GetErrorCode());
2858     return NULL;
2859   }
2860   aFillet->GetLastFunction()->SetDescription("");
2861
2862   TopoDS_Shape aFilletShape = aFillet->GetValue();
2863   aFunction->SetValue(aFilletShape);
2864   // END of fillet
2865
2866 // VSR: debug issues 0021568 and 0021550 (15/05/2012) - BEGIN (1)
2867 // the following block, when enabled, leads to partitioning problems
2868 #if 0
2869 // VSR: debug issues 0021568 and 0021550 (15/05/2012) - END (1)
2870   // BEGIN: Limit tolerances (debug)
2871   Handle(GEOM_Object) aCorr1 = myHealingOperations->LimitTolerance(aShape, 1e-07);
2872   TopoDS_Shape aCorr1Shape = aCorr1->GetValue();
2873   aShape->GetLastFunction()->SetValue(aCorr1Shape);
2874   aCorr1->GetLastFunction()->SetDescription("");
2875   // END: Limit tolerances (debug)
2876 // VSR: debug issues 0021568 and 0021550 (15/05/2012) - BEGIN (2)
2877 #endif
2878 // VSR: debug issues 0021568 and 0021550 (15/05/2012) - END (2)
2879
2880   if (theHexMesh) {
2881     if (!MakePipeTShapePartition(aShape, theR1, theW1, theL1, theR2, theW2, theL2, 0, 0, theRF, false))
2882       return NULL;
2883     if (!MakePipeTShapeMirrorAndGlue(aShape, theR1, theW1, theL1, theR2, theW2, theL2))
2884       return NULL;
2885   }
2886
2887   // Add thickness reduction elements
2888   // at the three extremities: Left, Right and Incident
2889   try {
2890 #if OCC_VERSION_LARGE > 0x06010000
2891     OCC_CATCH_SIGNALS;
2892 #endif
2893     if (isTRL || isTRR || isTRI) {
2894       TopoDS_Shape aResShape =
2895         MakePipeTShapeThicknessReduction(aShape->GetValue(), theR1, theW1, theL1, theR2, theW2, theL2,
2896                                          theRL, theWL, theLtransL, theLthinL,
2897                                          theRR, theWR, theLtransR, theLthinR,
2898                                          theRI, theWI, theLtransI, theLthinI,
2899                                          !theHexMesh);
2900       aFunction->SetValue(aResShape);
2901     }
2902   }
2903   catch (Standard_Failure) {
2904     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2905     SetErrorCode(aFail->GetMessageString());
2906     return NULL;
2907   }
2908
2909   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
2910   aSeq->Append(aShape);
2911
2912   try {
2913     if (theHexMesh) {
2914       // Get the groups
2915       if (!MakeGroups(aShape, TSHAPE_FILLET, theR1, theW1, theL1, theR2, theW2, theL2,
2916                       0., 0., theRF, aSeq, gp_Trsf()))
2917         return NULL;
2918     }
2919
2920     // Get internal group.
2921     if (!MakeInternalGroup(aShape, theR1, theL1, theR2, theL2, theRL, theLtransL,
2922                            theRR, theLtransR, theRI, theLtransI,
2923                            aSeq, gp_Trsf())) {
2924       return NULL;
2925     }
2926   }
2927   catch (Standard_Failure) {
2928     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
2929     SetErrorCode(aFail->GetMessageString());
2930     return NULL;
2931   }
2932
2933   //Make a Python command
2934   TCollection_AsciiString anEntry, aListRes("[");
2935   // Iterate over the sequence aSeq
2936   Standard_Integer aNbGroups = aSeq->Length();
2937   Standard_Integer i = 1;
2938   for (; i <= aNbGroups; i++) {
2939     Handle(Standard_Transient) anItem = aSeq->Value(i);
2940     if (anItem.IsNull()) continue;
2941     Handle(GEOM_Object) aGroup = Handle(GEOM_Object)::DownCast(anItem);
2942     if (aGroup.IsNull()) continue;
2943     //Make a Python command
2944     TDF_Tool::Entry(aGroup->GetEntry(), anEntry);
2945     aListRes += anEntry + ", ";
2946   }
2947   aListRes.Trunc(aListRes.Length() - 2);
2948
2949   GEOM::TPythonDump pd (aFunction);
2950
2951   pd << aListRes.ToCString() << "] = geompy.MakePipeTShapeFillet("
2952      << theR1 << ", " << theW1 << ", " << theL1 << ", "
2953      << theR2 << ", " << theW2 << ", " << theL2 << ", "
2954      << theRF << ", " << theHexMesh;
2955
2956   // thickness reduction
2957   if (isTRL)
2958     pd << ", theRL=" << theRL << ", theWL=" << theWL
2959        << ", theLtransL=" << theLtransL << ", theLthinL=" << theLthinL;
2960   if (isTRR)
2961     pd << ", theRR=" << theRR << ", theWR=" << theWR
2962        << ", theLtransR=" << theLtransR << ", theLthinR=" << theLthinR;
2963   if (isTRI)
2964     pd << ", theRI=" << theRI << ", theWI=" << theWI
2965        << ", theLtransI=" << theLtransI << ", theLthinI=" << theLthinI;
2966
2967   pd << ")";
2968
2969   SetErrorCode(OK);
2970
2971   return aSeq;
2972 }
2973
2974 //=============================================================================
2975 /*!
2976  *  MakePipeTShapeFilletWithPosition
2977  *  \brief Create a T-shape object with specified caracteristics for the main and
2978  *  the incident pipes (radius, width, half-length). A fillet is created
2979  *  on the junction of the pipes.
2980  *  The extremities of the main pipe are located on junctions points P1 and P2.
2981  *  The extremity of the incident pipe is located on junction point P3.
2982  *  \param theR1 Internal radius of main pipe
2983  *  \param theW1 Width of main pipe
2984  *  \param theL1 Half-length of main pipe
2985  *  \param theR2 Internal radius of incident pipe (R2 < R1)
2986  *  \param theW2 Width of incident pipe (R2+W2 < R1+W1)
2987  *  \param theL2 Half-length of incident pipe
2988  *  \param theRF Radius of curvature of fillet
2989  *  \param theHexMesh Boolean indicating if shape is prepared for hex mesh
2990  *  \param theP1 1st junction point of main pipe
2991  *  \param theP2 2nd junction point of main pipe
2992  *  \param theP3 Junction point of incident pipe
2993  *  \return List of GEOM_Objects, containing the created shape and propagation groups.
2994  */
2995 //=============================================================================
2996 Handle(TColStd_HSequenceOfTransient)
2997 GEOMImpl_IAdvancedOperations::MakePipeTShapeFilletWithPosition
2998                              (double theR1, double theW1, double theL1,
2999                               double theR2, double theW2, double theL2,
3000                               double theRL, double theWL, double theLtransL, double theLthinL,
3001                               double theRR, double theWR, double theLtransR, double theLthinR,
3002                               double theRI, double theWI, double theLtransI, double theLthinI,
3003                               double theRF, bool theHexMesh,
3004                               Handle(GEOM_Object) theP1,
3005                               Handle(GEOM_Object) theP2,
3006                               Handle(GEOM_Object) theP3)
3007 {
3008   SetErrorCode(KO);
3009   //Add a new object
3010   Handle(GEOM_Object) aShape = GetEngine()->AddObject(GetDocID(), GEOM_TSHAPE);
3011   //Add a new shape function with parameters
3012   Handle(GEOM_Function) aFunction = aShape->AddFunction(GEOMImpl_PipeTShapeDriver::GetID(), TSHAPE_FILLET);
3013   if (aFunction.IsNull()) return NULL;
3014
3015   //Check if the function is set correctly
3016   if (aFunction->GetDriverGUID() != GEOMImpl_PipeTShapeDriver::GetID()) return NULL;
3017
3018   // Check new position
3019   if (!CheckCompatiblePosition(theL1, theL2, theP1, theP2, theP3, 0.01)) {
3020     return NULL;
3021   }
3022
3023   GEOMImpl_IPipeTShape aData(aFunction);
3024
3025   aData.SetR1(theR1);
3026   aData.SetW1(theW1);
3027   aData.SetL1(theL1);
3028   aData.SetR2(theR2);
3029   aData.SetW2(theW2);
3030   aData.SetL2(theL2);
3031   aData.SetRF(theRF);
3032   aData.SetHexMesh(theHexMesh);
3033
3034   bool isTRL = (theRL + theWL + theLtransL + theLthinL) > Precision::Confusion();
3035   bool isTRR = (theRR + theWR + theLtransR + theLthinR) > Precision::Confusion();
3036   bool isTRI = (theRI + theWI + theLtransI + theLthinI) > Precision::Confusion();
3037
3038   //Compute the resulting value
3039   try {
3040 #if OCC_VERSION_LARGE > 0x06010000
3041     OCC_CATCH_SIGNALS;
3042 #endif
3043     if (!GetSolver()->ComputeFunction(aFunction)) {
3044       SetErrorCode("TShape driver failed");
3045       return NULL;
3046     }
3047   }
3048   catch (Standard_Failure) {
3049     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3050     SetErrorCode(aFail->GetMessageString());
3051     return NULL;
3052   }
3053
3054   // BEGIN of fillet
3055   TopoDS_Shape aShapeShape = aShape->GetValue();
3056   TopTools_IndexedMapOfShape anEdgesIndices;
3057   TopExp::MapShapes(aShapeShape, anEdgesIndices);
3058   // Common edges on external cylinders
3059   Handle(GEOM_Object) box_e;
3060   if (theHexMesh) {
3061     box_e = my3DPrimOperations->MakeBoxDXDYDZ(theR2+theW2, theR2+theW2, theR1+theW1);
3062   }
3063   else {
3064     box_e = my3DPrimOperations->MakeBoxDXDYDZ(2*(theR2+theW2), 2*(theR2+theW2), theR1+theW1);
3065   }
3066   box_e->GetLastFunction()->SetDescription("");
3067   box_e = myTransformOperations->TranslateDXDYDZ(box_e, -(theR2+theW2), -(theR2+theW2), 0);
3068   box_e->GetLastFunction()->SetDescription("");
3069
3070   Handle(TColStd_HSequenceOfInteger) edges_e =
3071     myShapesOperations->GetShapesOnBoxIDs(box_e, aShape, TopAbs_EDGE, GEOMAlgo_ST_IN);
3072   box_e->GetLastFunction()->SetDescription("");
3073
3074   if (edges_e.IsNull() || edges_e->Length() == 0) {
3075     SetErrorCode("External edges not found");
3076     return NULL;
3077   }
3078   int nbEdgesInFillet = 0;
3079   std::list<int> theEdges;
3080   for (int i=1; i<=edges_e->Length();i++) {
3081     int edgeID = edges_e->Value(i);
3082     TopoDS_Shape theEdge = anEdgesIndices.FindKey(edgeID);
3083     TopExp_Explorer Ex(theEdge,TopAbs_VERTEX);
3084     while (Ex.More()) {
3085       gp_Pnt aPt = BRep_Tool::Pnt(TopoDS::Vertex(Ex.Current()));
3086       if (Abs(aPt.Z() - (theR1+theW1)) <= Precision::Confusion()) {
3087         nbEdgesInFillet ++;
3088         theEdges.push_back(edgeID);
3089       }
3090       Ex.Next();
3091     }
3092     if (theHexMesh && nbEdgesInFillet == 1)
3093       break;
3094   }
3095
3096   Handle(GEOM_Object) aFillet;
3097   try {
3098     aFillet = myLocalOperations->MakeFilletEdges(aShape, theRF, theEdges);
3099   }
3100   catch (Standard_Failure) {
3101     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3102     SetErrorCode(aFail->GetMessageString());
3103     return NULL;
3104   }
3105   if (aFillet.IsNull()) {
3106     SetErrorCode("Fillet can not be computed on the given shape with the given parameters");
3107     return NULL;
3108   }
3109   aFillet->GetLastFunction()->SetDescription("");
3110
3111   TopoDS_Shape aFilletShape = aFillet->GetValue();
3112   aFunction->SetValue(aFilletShape);
3113   // END of fillet
3114
3115 // VSR: debug issues 0021568 and 0021550 (15/05/2012) - BEGIN (3)
3116 // the following block, when enabled, leads to partitioning problems
3117 #if 0
3118 // VSR: debug issues 0021568 and 0021550 (15/05/2012) - END (3)
3119   // BEGIN: Limit tolerances (debug)
3120   Handle(GEOM_Object) aCorr1 = myHealingOperations->LimitTolerance(aShape, 1e-07);
3121   TopoDS_Shape aCorr1Shape = aCorr1->GetValue();
3122   aShape->GetLastFunction()->SetValue(aCorr1Shape);
3123   aCorr1->GetLastFunction()->SetDescription("");
3124   // END: Limit tolerances (debug)
3125 // VSR: debug issues 0021568 and 0021550 (15/05/2012) - BEGIN (4)
3126 #endif
3127 // VSR: debug issues 0021568 and 0021550 (15/05/2012) - END (4)
3128
3129   if (theHexMesh) {
3130     if (!MakePipeTShapePartition(aShape, theR1, theW1, theL1, theR2, theW2, theL2, 0, 0, theRF, false))
3131       return NULL;
3132     if (!MakePipeTShapeMirrorAndGlue(aShape, theR1, theW1, theL1, theR2, theW2, theL2))
3133       return NULL;
3134   }
3135
3136   // Add thickness reduction elements
3137   // at the three extremities: Left, Right and Incident
3138   try {
3139 #if OCC_VERSION_LARGE > 0x06010000
3140     OCC_CATCH_SIGNALS;
3141 #endif
3142     if (isTRL || isTRR || isTRI) {
3143       TopoDS_Shape aResShape =
3144         MakePipeTShapeThicknessReduction(aShape->GetValue(), theR1, theW1, theL1, theR2, theW2, theL2,
3145                                          theRL, theWL, theLtransL, theLthinL,
3146                                          theRR, theWR, theLtransR, theLthinR,
3147                                          theRI, theWI, theLtransI, theLthinI,
3148                                          !theHexMesh);
3149       aFunction->SetValue(aResShape);
3150     }
3151   }
3152   catch (Standard_Failure) {
3153     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3154     SetErrorCode(aFail->GetMessageString());
3155     return NULL;
3156   }
3157
3158   // Set Position
3159   gp_Trsf aTrsf = GetPositionTrsf(theL1, theL2, theP1, theP2, theP3);
3160   BRepBuilderAPI_Transform aTransformation (aShape->GetValue(), aTrsf, Standard_False);
3161   TopoDS_Shape aTrsf_Shape = aTransformation.Shape();
3162   aFunction->SetValue(aTrsf_Shape);
3163
3164   Handle(TColStd_HSequenceOfTransient) aSeq = new TColStd_HSequenceOfTransient;
3165   aSeq->Append(aShape);
3166
3167   try {
3168     if (theHexMesh) {
3169       // Get the groups
3170       if (!MakeGroups(aShape, TSHAPE_FILLET, theR1, theW1, theL1, theR2, theW2, theL2,
3171                       0., 0., theRF, aSeq, aTrsf))
3172         return NULL;
3173     }
3174
3175     // Get internal group.
3176     if (!MakeInternalGroup(aShape, theR1, theL1, theR2, theL2, theRL, theLtransL,
3177                            theRR, theLtransR, theRI, theLtransI,
3178                            aSeq, aTrsf)) {
3179       return NULL;
3180     }
3181   }
3182   catch (Standard_Failure) {
3183     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3184     SetErrorCode(aFail->GetMessageString());
3185     return NULL;
3186   }
3187
3188   //Make a Python command
3189   TCollection_AsciiString anEntry, aListRes("[");
3190   // Iterate over the sequence aSeq
3191   Standard_Integer aNbGroups = aSeq->Length();
3192   Standard_Integer i = 1;
3193   for (; i <= aNbGroups; i++) {
3194     Handle(Standard_Transient) anItem = aSeq->Value(i);
3195     if (anItem.IsNull()) continue;
3196     Handle(GEOM_Object) aGroup = Handle(GEOM_Object)::DownCast(anItem);
3197     if (aGroup.IsNull()) continue;
3198     //Make a Python command
3199     TDF_Tool::Entry(aGroup->GetEntry(), anEntry);
3200     aListRes += anEntry + ", ";
3201   }
3202   aListRes.Trunc(aListRes.Length() - 2);
3203
3204   GEOM::TPythonDump pd (aFunction);
3205
3206   pd << aListRes.ToCString() << "] = geompy.MakePipeTShapeFillet("
3207      << theR1 << ", " << theW1 << ", " << theL1 << ", "
3208      << theR2  << ", " << theW2 << ", " << theL2 << ", "
3209      << theRF << ", " << theHexMesh << ", "
3210      << theP1 << ", " << theP2 << ", " << theP3;
3211
3212   // thickness reduction
3213   if (isTRL)
3214     pd << ", theRL=" << theRL << ", theWL=" << theWL
3215        << ", theLtransL=" << theLtransL << ", theLthinL=" << theLthinL;
3216   if (isTRR)
3217     pd << ", theRR=" << theRR << ", theWR=" << theWR
3218        << ", theLtransR=" << theLtransR << ", theLthinR=" << theLthinR;
3219   if (isTRI)
3220     pd << ", theRI=" << theRI << ", theWI=" << theWI
3221        << ", theLtransI=" << theLtransI << ", theLthinI=" << theLthinI;
3222
3223   pd << ")";
3224
3225   SetErrorCode(OK);
3226
3227   return aSeq;
3228 }
3229
3230 //=============================================================================
3231 /*!
3232  *  This function allows to create a disk already divided into blocks. It can be
3233  *  used to create divided pipes for later meshing in hexaedra.
3234  *  \param theR Radius of the disk
3235  *  \param theRatio Relative size of the central square diagonal against the disk diameter
3236  *  \param theOrientation Plane on which the disk will be built
3237  *  \param thePattern The division pattern of the disk (hexagon or square in the center)
3238  *  \return New GEOM_Object, containing the created shape.
3239  */
3240 //=============================================================================
3241 Handle(GEOM_Object) GEOMImpl_IAdvancedOperations::MakeDividedDisk (double theR, double theRatio,
3242                                                                    int theOrientation, int thePattern)
3243 {
3244   SetErrorCode(KO);
3245
3246   if (theOrientation != 1 &&
3247       theOrientation != 2 &&
3248       theOrientation != 3)
3249   {
3250     SetErrorCode("theOrientation must be 1(=OXY), 2(=OYZ) or 3(=OZX)");
3251     return NULL;
3252   }
3253   //Add a new object
3254   Handle(GEOM_Object) aShape = GetEngine()->AddObject(GetDocID(), GEOM_DIVIDEDDISK);
3255
3256   //Add a new shape function with parameters
3257   Handle(GEOM_Function) aFunction = aShape->AddFunction(GEOMImpl_DividedDiskDriver::GetID(), DIVIDEDDISK_R_RATIO);
3258   if (aFunction.IsNull()) return NULL;
3259
3260   //Check if the function is set correctly
3261   if (aFunction->GetDriverGUID() != GEOMImpl_DividedDiskDriver::GetID()) return NULL;
3262
3263   GEOMImpl_IDividedDisk aData (aFunction);
3264
3265   aData.SetR(theR);
3266   aData.SetRatio(theRatio);
3267   aData.SetOrientation(theOrientation);
3268   aData.SetType(thePattern);
3269
3270   //Compute the resulting value
3271   try {
3272 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
3273     OCC_CATCH_SIGNALS;
3274 #endif
3275     if (!GetSolver()->ComputeFunction(aFunction)) {
3276       SetErrorCode("DividedDisk driver failed");
3277       return NULL;
3278     }
3279   }
3280   catch (Standard_Failure) {
3281     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3282     SetErrorCode(aFail->GetMessageString());
3283     return NULL;
3284   }
3285
3286   std::string aPatternStr;
3287
3288   switch(thePattern)
3289   {
3290     case 0:
3291       aPatternStr = "GEOM.SQUARE";
3292       break;
3293     case 1:
3294       aPatternStr = "GEOM.HEXAGON";
3295       break;
3296   }
3297
3298   //Make a Python command
3299   GEOM::TPythonDump(aFunction) << aShape << " = geompy.MakeDividedDisk(" << theR << ", " << theOrientation << ", " << aPatternStr.c_str() << ")";
3300
3301   SetErrorCode(OK);
3302
3303   return aShape;
3304 }
3305
3306 //=============================================================================
3307 /*!
3308  *  This function allows to create a disk already divided into blocks. It can be
3309  *  used to create divided pipes for later meshing in hexaedra.
3310  *  \param theR Radius of the disk
3311  *  \param theRatio Relative size of the central square diagonal against the disk diameter
3312  *  \return New GEOM_Object, containing the created shape.
3313  */
3314 //=============================================================================
3315 Handle(GEOM_Object) GEOMImpl_IAdvancedOperations::MakeDividedDiskPntVecR (Handle(GEOM_Object) thePnt,
3316                                                                           Handle(GEOM_Object) theVec,
3317                                                                           double theR,
3318                                                                           double theRatio,
3319                                                                           int    thePattern)
3320 {
3321   SetErrorCode(KO);
3322
3323   //Add a new object
3324   Handle(GEOM_Object) aShape = GetEngine()->AddObject(GetDocID(), GEOM_DIVIDEDDISK);
3325
3326   //Add a new shape function with parameters
3327   Handle(GEOM_Function) aFunction = aShape->AddFunction(GEOMImpl_DividedDiskDriver::GetID(), DIVIDEDDISK_R_VECTOR_PNT);
3328   if (aFunction.IsNull()) return NULL;
3329
3330   //Check if the function is set correctly
3331   if (aFunction->GetDriverGUID() != GEOMImpl_DividedDiskDriver::GetID()) return NULL;
3332
3333   GEOMImpl_IDividedDisk aData (aFunction);
3334
3335   Handle(GEOM_Function) aRefPnt = thePnt->GetLastFunction();
3336   Handle(GEOM_Function) aRefVec = theVec->GetLastFunction();
3337
3338   if (aRefPnt.IsNull() || aRefVec.IsNull()) return NULL;
3339
3340   aData.SetCenter(aRefPnt);
3341   aData.SetVector(aRefVec);
3342
3343   aData.SetR(theR);
3344   aData.SetRatio(theRatio);
3345   aData.SetType(thePattern);
3346
3347   //Compute the resulting value
3348   try {
3349 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
3350     OCC_CATCH_SIGNALS;
3351 #endif
3352     if (!GetSolver()->ComputeFunction(aFunction)) {
3353       SetErrorCode("DividedDisk driver failed");
3354       return NULL;
3355     }
3356   }
3357   catch (Standard_Failure) {
3358     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3359     SetErrorCode(aFail->GetMessageString());
3360     return NULL;
3361   }
3362
3363   std::string aPatternStr;
3364
3365   switch(thePattern)
3366   {
3367     case 0:
3368       aPatternStr = "GEOM.SQUARE";
3369       break;
3370     case 1:
3371       aPatternStr = "GEOM.HEXAGON";
3372       break;
3373   }
3374
3375
3376   //Make a Python command
3377   GEOM::TPythonDump(aFunction) << aShape << " = geompy.MakeDividedDiskPntVecR(" << thePnt << ", " << theVec << ", " << theR << ", " << aPatternStr.c_str() << ")";
3378
3379   SetErrorCode(OK);
3380
3381   return aShape;
3382 }
3383
3384 //=============================================================================
3385 /*!
3386  *  Builds a cylinder prepared for hexa meshes
3387  *  \param theR Radius of the cylinder
3388  *  \param theH Height of the cylinder
3389  *  \return New GEOM_Object, containing the created shape.
3390  */
3391 //=============================================================================
3392 Handle(GEOM_Object) GEOMImpl_IAdvancedOperations::MakeDividedCylinder (double theR,
3393                                                                        double theH,
3394                                                                        int    thePattern)
3395 {
3396   SetErrorCode(KO);
3397
3398   //Add a new object
3399   Handle(GEOM_Object) aShape = GetEngine()->AddObject(GetDocID(), GEOM_DIVIDEDCYLINDER);
3400
3401   Handle(GEOM_Object) aBaseShape = MakeDividedDisk(theR, 67.0, 1, thePattern);
3402   aBaseShape->GetLastFunction()->SetDescription("");   // Erase dump of MakeDividedDisk
3403
3404   aShape = my3DPrimOperations->MakePrismDXDYDZ(aBaseShape,0.0,0.0,theH, -1.0);
3405
3406   Handle(GEOM_Function) aFunction =  aShape->GetLastFunction();
3407   aFunction->SetDescription("");   // Erase dump of MakePrismDXDYDZ
3408   aShape->SetType(GEOM_DIVIDEDCYLINDER);
3409
3410   std::string aPatternStr;
3411
3412   switch(thePattern)
3413   {
3414     case 0:
3415       aPatternStr = "GEOM.SQUARE";
3416       break;
3417     case 1:
3418       aPatternStr = "GEOM.HEXAGON";
3419       break;
3420   }
3421
3422   //Make a Python command
3423   GEOM::TPythonDump(aFunction) << aShape << " = geompy.MakeDividedCylinder(" << theR << ", " << theH << ", " << aPatternStr.c_str() << ")";
3424
3425   SetErrorCode(OK);
3426
3427   return aShape;
3428 }
3429 //=============================================================================
3430 /*!
3431  *  Create a smoothing surface from a set of points
3432  *  \param thelPoints list of points
3433  *  \return New GEOM_Object, containing the created shape.
3434  */
3435 //=============================================================================
3436 Handle(GEOM_Object) GEOMImpl_IAdvancedOperations::MakeSmoothingSurface (std::list<Handle(GEOM_Object)> thelPoints)
3437 {
3438   SetErrorCode(KO);
3439
3440   //Add a new object
3441   Handle(GEOM_Object) aShape = GetEngine()->AddObject(GetDocID(), GEOM_SMOOTHINGSURFACE);
3442
3443   //Add a new shape function with parameters
3444   Handle(GEOM_Function) aFunction = aShape->AddFunction(GEOMImpl_SmoothingSurfaceDriver::GetID(), SMOOTHINGSURFACE_LPOINTS);
3445   if (aFunction.IsNull()) return NULL;
3446
3447   //Check if the function is set correctly
3448   if (aFunction->GetDriverGUID() != GEOMImpl_SmoothingSurfaceDriver::GetID()) return NULL;
3449
3450   GEOMImpl_ISmoothingSurface aData (aFunction);
3451
3452   int aLen = thelPoints.size();
3453   aData.SetLength(aLen);
3454   int ind = 1;
3455   std::list<Handle(GEOM_Object)>::iterator it = thelPoints.begin();
3456   for (; it != thelPoints.end(); it++, ind++) {
3457     Handle(GEOM_Function) aRefPnt = (*it)->GetLastFunction();
3458     if (aRefPnt.IsNull()) {
3459       SetErrorCode("NULL point for bSplineFaceShape");
3460       return NULL;
3461     }
3462     aData.SetPoint(ind, aRefPnt);
3463   }
3464
3465
3466   //Compute the resulting value
3467   try {
3468 #if (OCC_VERSION_MAJOR << 16 | OCC_VERSION_MINOR << 8 | OCC_VERSION_MAINTENANCE) > 0x060100
3469     OCC_CATCH_SIGNALS;
3470 #endif
3471     if (!GetSolver()->ComputeFunction(aFunction)) {
3472       SetErrorCode("SmoothingSurface driver failed");
3473       return NULL;
3474     }
3475   }
3476   catch (Standard_Failure) {
3477     Handle(Standard_Failure) aFail = Standard_Failure::Caught();
3478     SetErrorCode(aFail->GetMessageString());
3479     return NULL;
3480   }
3481
3482   //Make a Python command
3483   GEOM::TPythonDump pd (aFunction);
3484   pd << aShape << " = geompy.MakeSmoothingSurface([";
3485   it = thelPoints.begin();
3486   pd << (*it++);
3487   while (it != thelPoints.end()) {
3488     pd << ", " << (*it++);
3489   }
3490   pd << "])";
3491
3492   SetErrorCode(OK);
3493
3494   return aShape;
3495 }
3496 /*@@ insert new functions before this line @@ do not remove this line @@ do not remove this line @@*/