src/GeomAPI/GeomAPI_Shell.cpp

   1 // Copyright (C) 2018-2023  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_Shell.h"
  21
  22 #include "GeomAPI_Ax3.h"
  23 #include "GeomAPI_Box.h"
  24 #include "GeomAPI_Cone.h"
  25 #include "GeomAPI_Cylinder.h"
  26 #include "GeomAPI_Face.h"
  27 #include "GeomAPI_Pnt.h"
  28 #include "GeomAPI_Sphere.h"
  29 #include "GeomAPI_Torus.h"
  30 #include "GeomAPI_Wire.h"
  31 #include "GeomAPI_XYZ.h"
  32
  33 #include <BRep_Builder.hxx>
  34 #include <BRepGProp.hxx>
  35 #include <GProp_GProps.hxx>
  36 #include <Precision.hxx>
  37 #include <TopExp_Explorer.hxx>
  38 #include <TopoDS.hxx>
  39
  40 #include <map>
  41
  42 //=================================================================================================
  43 GeomAPI_Shell::GeomAPI_Shell()
  44 {
  45   TopoDS_Shell* aShell = new TopoDS_Shell();
  46
  47   BRep_Builder aBuilder;
  48   aBuilder.MakeShell(*aShell);
  49
  50   this->setImpl(aShell);
  51 }
  52
  53 //=================================================================================================
  54 GeomAPI_Shell::GeomAPI_Shell(const std::shared_ptr<GeomAPI_Shape>& theShape)
  55 {
  56   if (!theShape->isNull() && theShape->isShell()) {
  57     setImpl(new TopoDS_Shape(theShape->impl<TopoDS_Shape>()));
  58   }
  59 }
  60
  61 //=================================================================================================
  62 std::shared_ptr<GeomAPI_Sphere> GeomAPI_Shell::getSphere() const
  63 {
  64   bool isSphere = true;
  65   bool isFirstFace = true;
  66   GeomSpherePtr aSphere;
  67
  68   for (TopExp_Explorer anExp(impl<TopoDS_Shape>(), TopAbs_FACE); anExp.More(); anExp.Next()) {
  69     GeomFacePtr aFace(new GeomAPI_Face);
  70     aFace->setImpl(new TopoDS_Shape(anExp.Current()));
  71
  72     GeomSpherePtr aCurSphere = aFace->getSphere();
  73     if (!aCurSphere) {
  74       isSphere = false;
  75       break;
  76     }
  77
  78     if (isFirstFace) {
  79       aSphere = aCurSphere;
  80       isFirstFace = false;
  81     }
  82     else if (aSphere->center()->distance(aCurSphere->center()) >= Precision::Confusion() ||
  83              Abs(aSphere->radius() - aCurSphere->radius()) >= Precision::Confusion()) {
  84       isSphere = false;
  85       break;
  86     }
  87   }
  88
  89   return isSphere ? aSphere : GeomSpherePtr();
  90 }
  91
  92 //=================================================================================================
  93 std::shared_ptr<GeomAPI_Cylinder> GeomAPI_Shell::getCylinder() const
  94 {
  95   bool isCylinder = true;
  96   bool isFirstFace = true;
  97
  98   GeomPointPtr aLocation;
  99   GeomDirPtr anAxis;
 100   double aRadius = 0.0;
 101   double aHeight = 0.0;
 102
 103   for (TopExp_Explorer anExp(impl<TopoDS_Shape>(), TopAbs_FACE); anExp.More(); anExp.Next()) {
 104     GeomFacePtr aFace(new GeomAPI_Face);
 105     aFace->setImpl(new TopoDS_Shape(anExp.Current()));
 106
 107     GeomCylinderPtr aCurCyl = aFace->getCylinder();
 108     if (!aCurCyl) {
 109       isCylinder = false;
 110       break;
 111     }
 112
 113     if (isFirstFace) {
 114       aLocation = aCurCyl->location();
 115       anAxis = aCurCyl->axis();
 116       aRadius = aCurCyl->radius();
 117       aHeight = aCurCyl->height();
 118       isFirstFace = false;
 119     }
 120     else {
 121       // compare radii
 122       if (Abs(aRadius - aCurCyl->radius()) >= Precision::Confusion() ||
 123         // check directions are collinear
 124         !anAxis->isParallel(aCurCyl->axis()) ||
 125         // check current center is on the main axis
 126         anAxis->xyz()->cross(aLocation->xyz()->decreased(aCurCyl->location()->xyz())
 127         )->squareModulus() >= Precision::SquareConfusion()) {
 128         isCylinder = false;
 129         break;
 130       }
 131
 132       double aMinHeight = 0.0;
 133       double aMaxHeight = aHeight;
 134
 135       std::shared_ptr<GeomAPI_XYZ> aCurCylLoc = aCurCyl->location()->xyz();
 136       std::shared_ptr<GeomAPI_XYZ> aCurCylLocHeight =
 137         aCurCylLoc->added(aCurCyl->axis()->xyz()->multiplied(aCurCyl->height()));
 138
 139       std::shared_ptr<GeomAPI_XYZ> anAxisXYZ = anAxis->xyz();
 140
 141       double aDist = anAxisXYZ->dot(aCurCylLoc->decreased(aLocation->xyz()));
 142       if (aDist < aMinHeight)
 143         aMinHeight = aDist;
 144       else if (aDist > aMaxHeight)
 145         aMaxHeight = aDist;
 146
 147       aDist = anAxisXYZ->dot(aCurCylLocHeight->decreased(aLocation->xyz()));
 148       if (aDist < aMinHeight)
 149         aMinHeight = aDist;
 150       else if (aDist > aMaxHeight)
 151         aMaxHeight = aDist;
 152
 153       if (aMinHeight < 0.0) {
 154         // move location of full cylinder
 155         aLocation->setX(aLocation->x() + aMinHeight * anAxis->x());
 156         aLocation->setY(aLocation->y() + aMinHeight * anAxis->y());
 157         aLocation->setZ(aLocation->z() + aMinHeight * anAxis->z());
 158       }
 159
 160       aHeight = aMaxHeight - aMinHeight;
 161     }
 162   }
 163
 164   GeomCylinderPtr aCylinder;
 165   if (isCylinder)
 166     aCylinder = GeomCylinderPtr(new GeomAPI_Cylinder(aLocation, anAxis, aRadius, aHeight));
 167   return aCylinder;
 168 }
 169
 170 //=================================================================================================
 171 std::shared_ptr<GeomAPI_Cone> GeomAPI_Shell::getCone() const
 172 {
 173   bool isCone = true;
 174   bool isFirstFace = true;
 175
 176   GeomPointPtr anApex;
 177   GeomDirPtr anAxis;
 178   double aSemiAngle = 0.0, aTanSemiAngle = 0.0;
 179   double aHeight1 = 0.0, aHeight2 = 0.0;
 180
 181   for (TopExp_Explorer anExp(impl<TopoDS_Shape>(), TopAbs_FACE); anExp.More(); anExp.Next()) {
 182     GeomFacePtr aFace(new GeomAPI_Face);
 183     aFace->setImpl(new TopoDS_Shape(anExp.Current()));
 184
 185     GeomConePtr aCurCone = aFace->getCone();
 186     if (!aCurCone) {
 187       isCone = false;
 188       break;
 189     }
 190
 191     if (isFirstFace) {
 192       anApex = aCurCone->apex();
 193       anAxis = aCurCone->axis();
 194       aSemiAngle = aCurCone->semiAngle();
 195       aTanSemiAngle = Tan(aSemiAngle);
 196       aHeight1 = aCurCone->radius1() / aTanSemiAngle;
 197       aHeight2 = aCurCone->radius2() / aTanSemiAngle;
 198       isFirstFace = false;
 199     }
 200     else {
 201       // check equal locations
 202       if (anApex->distance(aCurCone->apex()) >= Precision::Confusion() ||
 203       // check equal angles
 204           Abs(aSemiAngle - aCurCone->semiAngle()) >= Precision::Confusion() ||
 205       // check directions are collinear
 206           !anAxis->isParallel(aCurCone->axis())) {
 207         isCone = false;
 208         break;
 209       }
 210
 211       double aSign = anAxis->dot(aCurCone->axis());
 212       double aCurSemiAngle = aCurCone->semiAngle();
 213       double aTanCurSemiAngle = Tan(aCurSemiAngle);
 214
 215       double aH = aCurCone->radius1() / aTanCurSemiAngle * aSign;
 216       if (aH < aHeight1)
 217         aHeight1 = aH;
 218       else if (aH > aHeight2)
 219         aHeight2 = aH;
 220
 221       aH = aCurCone->radius2() / aTanCurSemiAngle * aSign;
 222       if (aH < aHeight1)
 223         aHeight1 = aH;
 224       else if (aH > aHeight2)
 225         aHeight2 = aH;
 226     }
 227   }
 228
 229   GeomConePtr aCone;
 230   if (isCone) {
 231     GeomPointPtr aLocation(new GeomAPI_Pnt(
 232       anApex->xyz()->added(anAxis->xyz()->multiplied(aHeight1))));
 233     double aRadius1 = aHeight1 * aTanSemiAngle;
 234     double aRadius2 = aHeight2 * aTanSemiAngle;
 235
 236     aCone = GeomConePtr(new GeomAPI_Cone(aLocation, anAxis, aSemiAngle, aRadius1, aRadius2));
 237   }
 238   return aCone;
 239 }
 240
 241 //=================================================================================================
 242 std::shared_ptr<GeomAPI_Torus> GeomAPI_Shell::getTorus() const
 243 {
 244   bool isTorus = true;
 245   bool isFirstFace = true;
 246   GeomTorusPtr aTorus;
 247
 248   for (TopExp_Explorer anExp(impl<TopoDS_Shape>(), TopAbs_FACE); anExp.More(); anExp.Next()) {
 249     GeomFacePtr aFace(new GeomAPI_Face);
 250     aFace->setImpl(new TopoDS_Shape(anExp.Current()));
 251
 252     GeomTorusPtr aCurTorus = aFace->getTorus();
 253     if (!aCurTorus) {
 254       isTorus = false;
 255       break;
 256     }
 257
 258     if (isFirstFace) {
 259       aTorus = aCurTorus;
 260       isFirstFace = false;
 261     }
 262     else {
 263       // compare radii
 264       if (Abs(aTorus->majorRadius() - aCurTorus->majorRadius()) >= Precision::Confusion() ||
 265           Abs(aTorus->minorRadius() - aCurTorus->minorRadius()) >= Precision::Confusion() ||
 266       // check equal centers
 267           aTorus->center()->distance(aCurTorus->center()) >= Precision::SquareConfusion() ||
 268       // check directions are collinear
 269           !aTorus->direction()->isParallel(aCurTorus->direction())) {
 270         isTorus = false;
 271         break;
 272       }
 273     }
 274   }
 275
 276   return isTorus ? aTorus : GeomTorusPtr();
 277 }
 278
 279 //=================================================================================================
 280 std::shared_ptr<GeomAPI_Box> GeomAPI_Shell::getParallelepiped() const
 281 {
 282   struct Plane
 283   {
 284     std::shared_ptr<GeomAPI_Ax3> myAxes;
 285     double myWidth;
 286     double myDepth;
 287     double myHeight;
 288   } aPlanes[6];
 289   std::map<int, int> aParallelPlanes;
 290
 291   int aNbPlanes = 0;
 292   for (TopExp_Explorer anExp(impl<TopoDS_Shape>(), TopAbs_WIRE);
 293        anExp.More() && aNbPlanes < 6; anExp.Next()) {
 294     GeomWirePtr aWire(new GeomAPI_Wire);
 295     aWire->setImpl(new TopoDS_Shape(anExp.Current()));
 296
 297     std::list<GeomPointPtr> aCorners;
 298     if (aWire->isRectangle(aCorners)) {
 299       // convert rectangle to plane with dimensions
 300
 301       // find corner with the smallest coordinates
 302       std::list<GeomPointPtr>::const_iterator aPrev = --aCorners.end();
 303       std::list<GeomPointPtr>::const_iterator aCur = aPrev--;
 304       std::list<GeomPointPtr>::const_iterator aNext = aCorners.begin();
 305       GeomPointPtr anOrigin = *aCur;
 306       GeomPointPtr aFront = *aNext;
 307       GeomPointPtr aBack = *aPrev;
 308       aPrev = aCur;
 309       aCur = aNext++;
 310       while (aNext != aCorners.end()) {
 311         if ((*aCur)->isLess(anOrigin)) {
 312           anOrigin = *aCur;
 313           aFront = *aNext;
 314           aBack = *aPrev;
 315         }
 316         aPrev = aCur;
 317         aCur = aNext++;
 318       }
 319
 320       aPlanes[aNbPlanes].myWidth = aBack->distance(anOrigin);
 321       aPlanes[aNbPlanes].myDepth = aFront->distance(anOrigin);
 322       aPlanes[aNbPlanes].myHeight = Precision::Infinite();
 323
 324       GeomDirPtr aDX(new GeomAPI_Dir(aBack->x() - anOrigin->x(),
 325                                      aBack->y() - anOrigin->y(),
 326                                      aBack->z() - anOrigin->z()));
 327       GeomDirPtr aDY(new GeomAPI_Dir(aFront->x() - anOrigin->x(),
 328                                      aFront->y() - anOrigin->y(),
 329                                      aFront->z() - anOrigin->z()));
 330       GeomDirPtr aDZ(new GeomAPI_Dir(aDX->cross(aDY)));
 331       aPlanes[aNbPlanes].myAxes =
 332           std::shared_ptr<GeomAPI_Ax3>(new GeomAPI_Ax3(anOrigin, aDX, aDZ));
 333
 334       // find parallel plane
 335       for (int i = 0; i < aNbPlanes; ++i) {
 336         double aDot = aPlanes[i].myAxes->normal()->dot(aDZ);
 337         if (Abs(aDot + 1.0) < Precision::Angular()) {
 338           if (aParallelPlanes.find(i) == aParallelPlanes.end())
 339             aParallelPlanes[i] = aNbPlanes;
 340           else
 341             break; // parallel planes already exist
 342         }
 343       }
 344
 345       ++aNbPlanes;
 346
 347     } else
 348       break;
 349   }
 350
 351   if (aNbPlanes != 6 || aParallelPlanes.size() != 3) // not a parallelepiped
 352     return GeomBoxPtr();
 353
 354   // calculate heights for planes computed by rectangles
 355   for (std::map<int, int>::iterator it = aParallelPlanes.begin();
 356        it != aParallelPlanes.end(); ++it) {
 357     GeomDirPtr aNormal = aPlanes[it->first].myAxes->normal();
 358     GeomPointPtr anOrigin = aPlanes[it->first].myAxes->origin();
 359     GeomPointPtr aNeighbor = aPlanes[it->second].myAxes->origin();
 360
 361     aPlanes[it->first].myHeight =
 362     aPlanes[it->second].myHeight =
 363         aNormal->xyz()->dot( aNeighbor->xyz()->decreased(anOrigin->xyz()) );
 364   }
 365
 366   // check if the box is oriented in the main axes
 367   int anIndex = 0;
 368   for (int i = 0; i < 6; ++i) {
 369     if (Abs(aPlanes[i].myAxes->dirX()->x() - 1.) < Precision::Angular() &&
 370         Abs(aPlanes[i].myAxes->normal()->z() - 1.) < Precision::Angular()) {
 371       anIndex = i;
 372       break;
 373     }
 374   }
 375
 376   // construct a box
 377   GeomBoxPtr aBox(new GeomAPI_Box(aPlanes[anIndex].myAxes,
 378                                   aPlanes[anIndex].myWidth,
 379                                   aPlanes[anIndex].myDepth,
 380                                   aPlanes[anIndex].myHeight));
 381   return aBox;
 382 }
 383
 384 //=================================================================================================
 385 GeomPointPtr GeomAPI_Shell::middlePoint() const
 386 {
 387   GeomPointPtr anInnerPoint;
 388
 389   const TopoDS_Shell& aShell = impl<TopoDS_Shell>();
 390   if (aShell.IsNull())
 391     return anInnerPoint;
 392
 393   GProp_GProps aProps;
 394   BRepGProp::SurfaceProperties(aShell, aProps, 1.e-4);
 395
 396   gp_Pnt aPnt = aProps.CentreOfMass();
 397   anInnerPoint = GeomPointPtr(new GeomAPI_Pnt(aPnt.X(), aPnt.Y(), aPnt.Z()));
 398   return anInnerPoint;
 399 }