src/GeomAPI/GeomAPI_Shell.cpp

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