src/GeomAPI/GeomAPI_Face.cpp

   1 // Copyright (C) 2014-2023  CEA, EDF
   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, or (at your option) any later version.
   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
  20 #include "GeomAPI_Face.h"
  21
  22 #include "GeomAPI_Dir.h"
  23 #include "GeomAPI_Pln.h"
  24 #include "GeomAPI_Pnt.h"
  25 #include "GeomAPI_Sphere.h"
  26 #include "GeomAPI_Curve.h"
  27 #include "GeomAPI_Cylinder.h"
  28 #include "GeomAPI_Cone.h"
  29 #include "GeomAPI_Torus.h"
  30
  31 #include <Basics_OCCTVersion.hxx>
  32
  33 #include <Bnd_Box2d.hxx>
  34 #include <BndLib_Add2dCurve.hxx>
  35 #include <BOPTools_AlgoTools.hxx>
  36 #include <BRep_Tool.hxx>
  37 #include <BRepAdaptor_Surface.hxx>
  38 #include <BRepGProp_Face.hxx>
  39 #include <BRepTools.hxx>
  40 #include <BRepTopAdaptor_TopolTool.hxx>
  41 #include <Geom_Surface.hxx>
  42 #include <Geom_SphericalSurface.hxx>
  43 #include <Geom_ConicalSurface.hxx>
  44 #include <Geom_CylindricalSurface.hxx>
  45 #include <Geom_OffsetSurface.hxx>
  46 #include <Geom_Plane.hxx>
  47 #include <Geom_RectangularTrimmedSurface.hxx>
  48 #include <Geom_SurfaceOfLinearExtrusion.hxx>
  49 #include <Geom_SurfaceOfRevolution.hxx>
  50 #include <Geom_SweptSurface.hxx>
  51 #include <Geom_ToroidalSurface.hxx>
  52 #include <GeomAPI_ExtremaCurveCurve.hxx>
  53 #include <GeomLib_IsPlanarSurface.hxx>
  54 #include <IntPatch_ImpImpIntersection.hxx>
  55 #include <IntTools_Context.hxx>
  56 #include <Standard_Type.hxx>
  57 #include <TopoDS.hxx>
  58 #include <TopoDS_Face.hxx>
  59
  60 #if OCC_VERSION_LARGE < 0x07070000
  61 #include <GeomAdaptor_HSurface.hxx>
  62 #else
  63 #include <GeomAdaptor_Surface.hxx>
  64 #endif
  65
  66 #include <gp_Sphere.hxx>
  67 #include <gp_Cylinder.hxx>
  68 #include <gp_Cone.hxx>
  69 #include <gp_Torus.hxx>
  70
  71 static void optimalBounds(const TopoDS_Face& theFace, double& theUMin, double& theUMax,
  72                                                       double& theVMin, double& theVMax);
  73
  74
  75 GeomAPI_Face::GeomAPI_Face()
  76   : GeomAPI_Shape()
  77 {
  78 }
  79
  80 GeomAPI_Face::GeomAPI_Face(const std::shared_ptr<GeomAPI_Shape>& theShape)
  81 {
  82   if (!theShape->isNull() && theShape->isFace()) {
  83     setImpl(new TopoDS_Shape(theShape->impl<TopoDS_Shape>()));
  84   }
  85 }
  86
  87 bool GeomAPI_Face::isEqual(std::shared_ptr<GeomAPI_Shape> theFace) const
  88 {
  89   if (!theFace.get())
  90     return false;
  91
  92   if (!theFace->isFace())
  93     return false;
  94
  95   const TopoDS_Shape& aMyShape = const_cast<GeomAPI_Face*>(this)->impl<TopoDS_Shape>();
  96   const TopoDS_Shape& aInShape = theFace->impl<TopoDS_Shape>();
  97
  98   TopoDS_Face aMyFace = TopoDS::Face(aMyShape);
  99   TopoDS_Face aInFace = TopoDS::Face(aInShape);
 100
 101   Handle(Geom_Surface) aMySurf = BRep_Tool::Surface(aMyFace);
 102   Handle(Geom_Surface) aInSurf = BRep_Tool::Surface(aInFace);
 103
 104   // Check that surfaces a the same type
 105   if (aMySurf->DynamicType() != aInSurf->DynamicType())
 106     return false;
 107
 108   // Get parameters of surfaces
 109   double aMyUMin, aMyUMax, aMyVMin, aMyVMax;
 110   aMySurf->Bounds(aMyUMin, aMyUMax, aMyVMin, aMyVMax);
 111   double aInUMin, aInUMax, aInVMin, aInVMax;
 112   aInSurf->Bounds(aInUMin, aInUMax, aInVMin, aInVMax);
 113
 114   // Check that parameters are the same
 115   if (fabs(aMyUMin - aInUMin) > Precision::PConfusion() ||
 116       fabs(aMyUMax - aInUMax) > Precision::PConfusion() ||
 117       fabs(aMyVMin - aInVMin) > Precision::PConfusion() ||
 118       fabs(aMyVMax - aInVMax) > Precision::PConfusion())
 119     return false;
 120
 121   Handle(IntTools_Context) aContext = new IntTools_Context();
 122   // Double check needed because BOPTools_AlgoTools::AreFacesSameDomain not very smart.
 123   Standard_Boolean aRes = BOPTools_AlgoTools::AreFacesSameDomain(aMyFace, aInFace, aContext)
 124     && BOPTools_AlgoTools::AreFacesSameDomain(aInFace, aMyFace, aContext);
 125
 126   return aRes == Standard_True;
 127 }
 128
 129 static Handle(Geom_Surface) baseSurface(const TopoDS_Face& theFace)
 130 {
 131   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(theFace);
 132   while (aSurf->IsKind(STANDARD_TYPE(Geom_RectangularTrimmedSurface))) {
 133     Handle(Geom_RectangularTrimmedSurface) rts =
 134         Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurf);
 135     aSurf = rts->BasisSurface();
 136   }
 137   return aSurf;
 138 }
 139
 140 bool GeomAPI_Face::isSameGeometry(const std::shared_ptr<GeomAPI_Shape> theShape) const
 141 {
 142   if (!theShape->isFace())
 143     return false;
 144   if (isSame(theShape))
 145     return true;
 146
 147   GeomFacePtr anOther = theShape->face();
 148   if (isPlanar() && anOther->isPlanar()) {
 149     GeomPlanePtr anOwnPlane = getPlane();
 150     GeomPlanePtr anOtherPlane = anOther->getPlane();
 151     return anOwnPlane->isCoincident(anOtherPlane);
 152   }
 153
 154   TopoDS_Face anOwnFace = TopoDS::Face(impl<TopoDS_Shape>());
 155   TopoDS_Face anOtherFace = TopoDS::Face(theShape->impl<TopoDS_Shape>());
 156
 157   Handle(Geom_Surface) anOwnSurf = baseSurface(anOwnFace);
 158   Handle(Geom_Surface) anOtherSurf = baseSurface(anOtherFace);
 159   if (anOwnSurf == anOtherSurf)
 160     return true;
 161
 162   // case of two elementary surfaces
 163   if (anOwnSurf->IsKind(STANDARD_TYPE(Geom_ElementarySurface)) &&
 164       anOtherSurf->IsKind(STANDARD_TYPE(Geom_ElementarySurface)))
 165   {
 166 #if OCC_VERSION_LARGE < 0x07070000
 167     Handle(GeomAdaptor_HSurface) aGA1 = new GeomAdaptor_HSurface(anOwnSurf);
 168     Handle(GeomAdaptor_HSurface) aGA2 = new GeomAdaptor_HSurface(anOtherSurf);
 169 #else
 170     Handle(GeomAdaptor_Surface) aGA1 = new GeomAdaptor_Surface(anOwnSurf);
 171     Handle(GeomAdaptor_Surface) aGA2 = new GeomAdaptor_Surface(anOtherSurf);
 172 #endif
 173
 174     Handle(BRepTopAdaptor_TopolTool) aTT1 = new BRepTopAdaptor_TopolTool();
 175     Handle(BRepTopAdaptor_TopolTool) aTT2 = new BRepTopAdaptor_TopolTool();
 176
 177     try {
 178       IntPatch_ImpImpIntersection anIIInt(aGA1, aTT1, aGA2, aTT2,
 179                                           Precision::Confusion(),
 180                                           Precision::Confusion());
 181       if (!anIIInt.IsDone() || anIIInt.IsEmpty())
 182         return false;
 183
 184       return anIIInt.TangentFaces();
 185     }
 186     catch (Standard_Failure const&) {
 187       return false;
 188     }
 189   }
 190
 191   // case of two cylindrical surfaces, at least one of which is a swept surface
 192   // swept surfaces: SurfaceOfLinearExtrusion, SurfaceOfRevolution
 193   if ((anOwnSurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)) ||
 194        anOwnSurf->IsKind(STANDARD_TYPE(Geom_SweptSurface))) &&
 195       (anOtherSurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)) ||
 196        anOtherSurf->IsKind(STANDARD_TYPE(Geom_SweptSurface))))
 197   {
 198     GeomCylinderPtr anOwnCyl = getCylinder();
 199     GeomCylinderPtr anOtherCyl = anOther->getCylinder();
 200     if (anOwnCyl && anOtherCyl)
 201       return anOwnCyl->isCoincident(anOtherCyl);
 202
 203     // compare two swept surfaces of the same type
 204     if ((anOwnSurf->IsKind(STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion)) &&
 205          anOtherSurf->IsKind(STANDARD_TYPE(Geom_SurfaceOfLinearExtrusion))) ||
 206         (anOwnSurf->IsKind(STANDARD_TYPE(Geom_SurfaceOfRevolution)) &&
 207          anOtherSurf->IsKind(STANDARD_TYPE(Geom_SurfaceOfRevolution)))) {
 208       Handle(Geom_SweptSurface) anOwnSwept = Handle(Geom_SweptSurface)::DownCast(anOwnSurf);
 209       Handle(Geom_SweptSurface) anOtherSwept = Handle(Geom_SweptSurface)::DownCast(anOtherSurf);
 210
 211       const gp_Dir& anOwnDir = anOwnSwept->Direction();
 212       const gp_Dir& anOtherDir = anOtherSwept->Direction();
 213
 214       if (anOwnDir.IsParallel(anOtherDir, Precision::Angular())) {
 215         Handle(Geom_Curve) anOwnCurve = anOwnSwept->BasisCurve();
 216         Handle(Geom_Curve) anOtherCurve = anOtherSwept->BasisCurve();
 217         GeomAPI_ExtremaCurveCurve anExtrema(anOwnCurve, anOtherCurve);
 218         return anExtrema.Extrema().IsParallel() &&
 219                anExtrema.TotalLowerDistance() < Precision::Confusion();
 220       }
 221     }
 222   }
 223
 224   return false;
 225 }
 226
 227 bool GeomAPI_Face::isCylindrical() const
 228 {
 229   const TopoDS_Shape& aShape = const_cast<GeomAPI_Face*>(this)->impl<TopoDS_Shape>();
 230   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(TopoDS::Face(aShape));
 231   Handle(Geom_RectangularTrimmedSurface) aTrimmed =
 232     Handle(Geom_RectangularTrimmedSurface)::DownCast(aSurf);
 233   if (!aTrimmed.IsNull())
 234     aSurf = aTrimmed->BasisSurface();
 235   return aSurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface)) == Standard_True;
 236 }
 237
 238 std::shared_ptr<GeomAPI_Pln> GeomAPI_Face::getPlane() const
 239 {
 240   std::shared_ptr<GeomAPI_Pln> aResult;
 241   TopoDS_Shape aShape = this->impl<TopoDS_Shape>();
 242   if (aShape.IsNull())
 243     return aResult;  // null shape
 244   if (aShape.ShapeType() != TopAbs_FACE)
 245     return aResult;  // not face
 246   TopoDS_Face aFace = TopoDS::Face(aShape);
 247   if (aFace.IsNull())
 248     return aResult;  // not face
 249   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
 250   if (aSurf.IsNull())
 251     return aResult;  // no surface
 252   GeomLib_IsPlanarSurface isPlanarSurf(aSurf);
 253   gp_Pln aPln;
 254   bool isPlanar = false;
 255   if (isPlanarSurf.IsPlanar()) {
 256     aPln = isPlanarSurf.Plan();
 257     isPlanar = true;
 258   }
 259   else if (aSurf->IsKind(STANDARD_TYPE(Geom_OffsetSurface))) {
 260     Handle(Geom_OffsetSurface) anOffsetSurf = Handle(Geom_OffsetSurface)::DownCast(aSurf);
 261     Handle(Geom_Surface) aBasisSurf = anOffsetSurf->BasisSurface();
 262     if (aBasisSurf->IsKind(STANDARD_TYPE(Geom_Plane))) {
 263       aPln = Handle(Geom_Plane)::DownCast(aBasisSurf)->Pln();
 264       gp_Vec aTranslation(aPln.Axis().Direction().XYZ() * anOffsetSurf->Offset());
 265       aPln.Translate(aTranslation);
 266       isPlanar = true;
 267     }
 268   }
 269
 270   if (isPlanar) {
 271     double aA, aB, aC, aD;
 272     aPln.Coefficients(aA, aB, aC, aD);
 273     if (aFace.Orientation() == TopAbs_REVERSED) {
 274       aA = -aA;
 275       aB = -aB;
 276       aC = -aC;
 277       aD = -aD;
 278     }
 279     aResult = std::shared_ptr<GeomAPI_Pln>(new GeomAPI_Pln(aA, aB, aC, aD));
 280   }
 281   return aResult;
 282 }
 283
 284 std::shared_ptr<GeomAPI_Sphere> GeomAPI_Face::getSphere() const
 285 {
 286   GeomSpherePtr aSphere;
 287
 288   const TopoDS_Face& aFace = TopoDS::Face(impl<TopoDS_Shape>());
 289   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
 290   if (aSurf->IsKind(STANDARD_TYPE(Geom_SphericalSurface))) {
 291     gp_Sphere aSph = Handle(Geom_SphericalSurface)::DownCast(aSurf)->Sphere();
 292     const gp_Pnt& aCenter = aSph.Location();
 293     double aRadius = aSph.Radius();
 294
 295     GeomPointPtr aCenterPnt(new GeomAPI_Pnt(aCenter.X(), aCenter.Y(), aCenter.Z()));
 296     aSphere = GeomSpherePtr(new GeomAPI_Sphere(aCenterPnt, aRadius));
 297   }
 298   return aSphere;
 299 }
 300
 301 std::shared_ptr<GeomAPI_Cylinder> GeomAPI_Face::getCylinder() const
 302 {
 303   GeomCylinderPtr aCylinder;
 304
 305   const TopoDS_Face& aFace = TopoDS::Face(impl<TopoDS_Shape>());
 306   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
 307   if (aSurf->IsKind(STANDARD_TYPE(Geom_CylindricalSurface))) {
 308     gp_Cylinder aCyl = Handle(Geom_CylindricalSurface)::DownCast(aSurf)->Cylinder();
 309     gp_Pnt aLoc = aCyl.Location();
 310     const gp_Dir& aDir = aCyl.Position().Direction();
 311     double aRadius = aCyl.Radius();
 312
 313     double aUMin, aUMax, aVMin, aVMax;
 314     BRepTools::UVBounds(aFace, aUMin, aUMax, aVMin, aVMax);
 315     double aHeight = aVMax - aVMin;
 316
 317     aLoc.ChangeCoord() += aDir.XYZ() * aVMin;
 318     GeomPointPtr aLocation(new GeomAPI_Pnt(aLoc.X(), aLoc.Y(), aLoc.Z()));
 319     GeomDirPtr aDirection(new GeomAPI_Dir(aDir.X(), aDir.Y(), aDir.Z()));
 320     aCylinder = GeomCylinderPtr(new GeomAPI_Cylinder(aLocation, aDirection, aRadius, aHeight));
 321   }
 322   return aCylinder;
 323 }
 324
 325 std::shared_ptr<GeomAPI_Cone> GeomAPI_Face::getCone() const
 326 {
 327   GeomConePtr aCone;
 328
 329   const TopoDS_Face& aFace = TopoDS::Face(impl<TopoDS_Shape>());
 330   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
 331   if (aSurf->IsKind(STANDARD_TYPE(Geom_ConicalSurface))) {
 332     gp_Cone aCon = Handle(Geom_ConicalSurface)::DownCast(aSurf)->Cone();
 333     gp_Pnt aLoc = aCon.Location();
 334     gp_Dir aDir = aCon.Position().Direction();
 335
 336     double aUMin, aUMax, aVMin, aVMax;
 337     BRepTools::UVBounds(aFace, aUMin, aUMax, aVMin, aVMax);
 338
 339     double aSemiAngle = Abs(aCon.SemiAngle());
 340     double aRadius1 = Abs(aCon.RefRadius() + aVMin * Sin(aCon.SemiAngle()));
 341     double aRadius2 = Abs(aCon.RefRadius() + aVMax * Sin(aCon.SemiAngle()));
 342
 343     aLoc.ChangeCoord() += aDir.XYZ() * aVMin * Cos(aCon.SemiAngle());
 344
 345     GeomPointPtr aLocation(new GeomAPI_Pnt(aLoc.X(), aLoc.Y(), aLoc.Z()));
 346     GeomDirPtr aDirection(new GeomAPI_Dir(aDir.X(), aDir.Y(), aDir.Z()));
 347     aCone = GeomConePtr(new GeomAPI_Cone(aLocation, aDirection, aSemiAngle, aRadius1, aRadius2));
 348   }
 349   return aCone;
 350 }
 351
 352 std::shared_ptr<GeomAPI_Torus> GeomAPI_Face::getTorus() const
 353 {
 354   GeomTorusPtr aTorus;
 355
 356   const TopoDS_Face& aFace = TopoDS::Face(impl<TopoDS_Shape>());
 357   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
 358   if (aSurf->IsKind(STANDARD_TYPE(Geom_ToroidalSurface))) {
 359     gp_Torus aTor = Handle(Geom_ToroidalSurface)::DownCast(aSurf)->Torus();
 360     const gp_Pnt& aLoc = aTor.Location();
 361     const gp_Dir& aDir = aTor.Position().Direction();
 362     double aMajorRadius = aTor.MajorRadius();
 363     double aMinorRadius = aTor.MinorRadius();
 364
 365     GeomPointPtr aCenter(new GeomAPI_Pnt(aLoc.X(), aLoc.Y(), aLoc.Z()));
 366     GeomDirPtr aDirection(new GeomAPI_Dir(aDir.X(), aDir.Y(), aDir.Z()));
 367     aTorus = GeomTorusPtr(new GeomAPI_Torus(aCenter, aDirection, aMajorRadius, aMinorRadius));
 368   }
 369   return aTorus;
 370 }
 371
 372 GeomPointPtr GeomAPI_Face::middlePoint() const
 373 {
 374   GeomPointPtr anInnerPoint;
 375
 376   const TopoDS_Face& aFace = impl<TopoDS_Face>();
 377   if (aFace.IsNull())
 378     return anInnerPoint;
 379
 380   double aUMin, aUMax, aVMin, aVMax;
 381   optimalBounds(aFace, aUMin, aUMax, aVMin, aVMax);
 382
 383   Handle(Geom_Surface) aSurf = BRep_Tool::Surface(aFace);
 384   if (aSurf.IsNull())
 385     return anInnerPoint;
 386
 387   gp_Pnt aPnt = aSurf->Value((aUMin + aUMax) * 0.5, (aVMin + aVMax) * 0.5);
 388   anInnerPoint = GeomPointPtr(new GeomAPI_Pnt(aPnt.X(), aPnt.Y(), aPnt.Z()));
 389   return anInnerPoint;
 390 }
 391
 392
 393 // ==================     Auxiliary functions     ========================
 394
 395 void optimalBounds(const TopoDS_Face& theFace, const TopoDS_Edge& theEdge, Bnd_Box2d& theBndBox)
 396 {
 397   Standard_Real aFirst, aLast;
 398   const Handle(Geom2d_Curve) aC2D = BRep_Tool::CurveOnSurface(theEdge, theFace, aFirst, aLast);
 399   if (aC2D.IsNull())
 400     return;
 401
 402   Standard_Real aXmin = 0.0, aYmin = 0.0, aXmax = 0.0, aYmax = 0.0;
 403   Standard_Real aUmin, aUmax, aVmin, aVmax;
 404   Bnd_Box2d aBoxC, aBoxS;
 405   BndLib_Add2dCurve::AddOptimal(aC2D, aFirst, aLast, 0., aBoxC);
 406   if (aBoxC.IsVoid())
 407     return;
 408
 409   aBoxC.Get(aXmin, aYmin, aXmax, aYmax);
 410
 411   TopLoc_Location aLoc;
 412   Handle(Geom_Surface) aS = BRep_Tool::Surface(theFace, aLoc);
 413   aS->Bounds(aUmin, aUmax, aVmin, aVmax);
 414
 415   if (aS->DynamicType() == STANDARD_TYPE(Geom_RectangularTrimmedSurface))
 416   {
 417     const Handle(Geom_RectangularTrimmedSurface) aSt =
 418         Handle(Geom_RectangularTrimmedSurface)::DownCast(aS);
 419     aS = aSt->BasisSurface();
 420   }
 421
 422   //
 423   if (!aS->IsUPeriodic())
 424   {
 425     Standard_Boolean isUPeriodic = Standard_False;
 426
 427     // Additional verification for U-periodicity for B-spline surfaces.
 428     // 1. Verify that the surface is U-closed (if such flag is false). Verification uses 2 points.
 429     // 2. Verify periodicity of surface inside UV-bounds of the edge. It uses 3 or 6 points.
 430     if (aS->DynamicType() == STANDARD_TYPE(Geom_BSplineSurface) &&
 431       (aXmin < aUmin || aXmax > aUmax))
 432     {
 433       Standard_Real aTol2 = 100 * Precision::SquareConfusion();
 434       isUPeriodic = Standard_True;
 435       gp_Pnt P1, P2;
 436       // 1. Verify that the surface is U-closed
 437       if (!aS->IsUClosed())
 438       {
 439         Standard_Real aVStep = aVmax - aVmin;
 440         for (Standard_Real aV = aVmin; aV <= aVmax; aV += aVStep)
 441         {
 442           P1 = aS->Value(aUmin, aV);
 443           P2 = aS->Value(aUmax, aV);
 444           if (P1.SquareDistance(P2) > aTol2)
 445           {
 446             isUPeriodic = Standard_False;
 447             break;
 448           }
 449         }
 450       }
 451       // 2. Verify periodicity of surface inside UV-bounds of the edge
 452       if (isUPeriodic) // the flag still not changed
 453       {
 454         Standard_Real aV = (aVmin + aVmax) * 0.5;
 455         Standard_Real aU[6]; // values of U lying out of surface boundaries
 456         Standard_Real aUpp[6]; // corresponding U-values plus/minus period
 457         Standard_Integer aNbPnt = 0;
 458         if (aXmin < aUmin)
 459         {
 460           aU[0] = aXmin;
 461           aU[1] = (aXmin + aUmin) * 0.5;
 462           aU[2] = aUmin;
 463           aUpp[0] = aU[0] + aUmax - aUmin;
 464           aUpp[1] = aU[1] + aUmax - aUmin;
 465           aUpp[2] = aU[2] + aUmax - aUmin;
 466           aNbPnt += 3;
 467         }
 468         if (aXmax > aUmax)
 469         {
 470           aU[aNbPnt] = aUmax;
 471           aU[aNbPnt + 1] = (aXmax + aUmax) * 0.5;
 472           aU[aNbPnt + 2] = aXmax;
 473           aUpp[aNbPnt] = aU[aNbPnt] - aUmax + aUmin;
 474           aUpp[aNbPnt + 1] = aU[aNbPnt + 1] - aUmax + aUmin;
 475           aUpp[aNbPnt + 2] = aU[aNbPnt + 2] - aUmax + aUmin;
 476           aNbPnt += 3;
 477         }
 478         for (Standard_Integer anInd = 0; anInd < aNbPnt; anInd++)
 479         {
 480           P1 = aS->Value(aU[anInd], aV);
 481           P2 = aS->Value(aUpp[anInd], aV);
 482           if (P1.SquareDistance(P2) > aTol2)
 483           {
 484             isUPeriodic = Standard_False;
 485             break;
 486           }
 487         }
 488       }
 489     }
 490
 491     if (!isUPeriodic)
 492     {
 493       if ((aXmin < aUmin) && (aUmin < aXmax))
 494       {
 495         aXmin = aUmin;
 496       }
 497       if ((aXmin < aUmax) && (aUmax < aXmax))
 498       {
 499         aXmax = aUmax;
 500       }
 501     }
 502   }
 503
 504   if (!aS->IsVPeriodic())
 505   {
 506     Standard_Boolean isVPeriodic = Standard_False;
 507
 508     // Additional verification for V-periodicity for B-spline surfaces.
 509     // 1. Verify that the surface is V-closed (if such flag is false). Verification uses 2 points.
 510     // 2. Verify periodicity of surface inside UV-bounds of the edge. It uses 3 or 6 points.
 511     if (aS->DynamicType() == STANDARD_TYPE(Geom_BSplineSurface) &&
 512       (aYmin < aVmin || aYmax > aVmax))
 513     {
 514       Standard_Real aTol2 = 100 * Precision::SquareConfusion();
 515       isVPeriodic = Standard_True;
 516       gp_Pnt P1, P2;
 517       // 1. Verify that the surface is V-closed
 518       if (!aS->IsVClosed())
 519       {
 520         Standard_Real aUStep = aUmax - aUmin;
 521         for (Standard_Real aU = aUmin; aU <= aUmax; aU += aUStep)
 522         {
 523           P1 = aS->Value(aU, aVmin);
 524           P2 = aS->Value(aU, aVmax);
 525           if (P1.SquareDistance(P2) > aTol2)
 526           {
 527             isVPeriodic = Standard_False;
 528             break;
 529           }
 530         }
 531       }
 532       // 2. Verify periodicity of surface inside UV-bounds of the edge
 533       if (isVPeriodic) // the flag still not changed
 534       {
 535         Standard_Real aU = (aUmin + aUmax) * 0.5;
 536         Standard_Real aV[6]; // values of V lying out of surface boundaries
 537         Standard_Real aVpp[6]; // corresponding V-values plus/minus period
 538         Standard_Integer aNbPnt = 0;
 539         if (aYmin < aVmin)
 540         {
 541           aV[0] = aYmin;
 542           aV[1] = (aYmin + aVmin) * 0.5;
 543           aV[2] = aVmin;
 544           aVpp[0] = aV[0] + aVmax - aVmin;
 545           aVpp[1] = aV[1] + aVmax - aVmin;
 546           aVpp[2] = aV[2] + aVmax - aVmin;
 547           aNbPnt += 3;
 548         }
 549         if (aYmax > aVmax)
 550         {
 551           aV[aNbPnt] = aVmax;
 552           aV[aNbPnt + 1] = (aYmax + aVmax) * 0.5;
 553           aV[aNbPnt + 2] = aYmax;
 554           aVpp[aNbPnt] = aV[aNbPnt] - aVmax + aVmin;
 555           aVpp[aNbPnt + 1] = aV[aNbPnt + 1] - aVmax + aVmin;
 556           aVpp[aNbPnt + 2] = aV[aNbPnt + 2] - aVmax + aVmin;
 557           aNbPnt += 3;
 558         }
 559         for (Standard_Integer anInd = 0; anInd < aNbPnt; anInd++)
 560         {
 561           P1 = aS->Value(aU, aV[anInd]);
 562           P2 = aS->Value(aU, aVpp[anInd]);
 563           if (P1.SquareDistance(P2) > aTol2)
 564           {
 565             isVPeriodic = Standard_False;
 566             break;
 567           }
 568         }
 569       }
 570     }
 571
 572     if (!isVPeriodic)
 573     {
 574       if ((aYmin < aVmin) && (aVmin < aYmax))
 575       {
 576         aYmin = aVmin;
 577       }
 578       if ((aYmin < aVmax) && (aVmax < aYmax))
 579       {
 580         aYmax = aVmax;
 581       }
 582     }
 583   }
 584
 585   aBoxS.Update(aXmin, aYmin, aXmax, aYmax);
 586
 587   theBndBox.Add(aBoxS);
 588 }
 589
 590 void optimalBounds(const TopoDS_Face& theFace, double& theUMin, double& theUMax,
 591                                                double& theVMin, double& theVMax)
 592 {
 593   Bnd_Box2d aBB;
 594
 595   for (TopExp_Explorer anExp(theFace, TopAbs_EDGE); anExp.More(); anExp.Next())
 596     optimalBounds(theFace, TopoDS::Edge(anExp.Current()), aBB);
 597
 598   aBB.Get(theUMin, theVMin, theUMax, theVMax);
 599 }