src/GeomAPI/GeomAPI_Face.cpp

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