G4OCCTSolid.hh

Go to the documentation of this file.

// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2026 G4OCCT Contributors
 
 
#ifndef G4OCCT_G4OCCTSolid_hh
#define G4OCCT_G4OCCTSolid_hh
 
#include <G4Cache.hh>
#include <G4Polyhedron.hh>
#include <G4ThreeVector.hh>
#include <G4VSolid.hh>
 
// OCCT shape representation
#include <BRepAdaptor_Surface.hxx>
#include <BRepClass3d_SolidClassifier.hxx>
#include <BRepExtrema_TriangleSet.hxx>
#include <Bnd_Box.hxx>
#include <IntCurvesFace_Intersector.hxx>
#include <TopoDS_Face.hxx>
#include <TopoDS_Shape.hxx>
#include <TopoDS_TShape.hxx>
#include <gp_Pln.hxx>
#include <gp_Pnt2d.hxx>
 
#include <atomic>
#include <condition_variable>
#include <cstdint>
#include <limits>
#include <memory>
#include <mutex>
#include <optional>
#include <stdexcept>
#include <string>
#include <unordered_map>
#include <utility>
#include <vector>
 
class G4OCCTSolid : public G4VSolid {
public:
  G4OCCTSolid(const G4String& name, const TopoDS_Shape& shape);
 
  ~G4OCCTSolid() override;
 
  static G4OCCTSolid* FromSTEP(const G4String& name, const std::string& path);
 
  // ── G4VSolid pure-virtual interface ───────────────────────────────────────
 
  EInside Inside(const G4ThreeVector& p) const override;
 
  G4ThreeVector SurfaceNormal(const G4ThreeVector& p) const override;
 
  G4double DistanceToIn(const G4ThreeVector& p, const G4ThreeVector& v) const override;
 
  G4double DistanceToIn(const G4ThreeVector& p) const override;
 
  G4double DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
                         const G4bool calcNorm = false, G4bool* validNorm = nullptr,
                         G4ThreeVector* n = nullptr) const override;
 
  G4double DistanceToOut(const G4ThreeVector& p) const override;
 
  G4ThreeVector GetPointOnSurface() const override;
 
  // ── G4OCCTSolid distance functions ────────────────────────────────────────
 
  G4double ExactDistanceToIn(const G4ThreeVector& p) const;
 
  G4double ExactDistanceToOut(const G4ThreeVector& p) const;
 
  G4double GetCubicVolume() override;
 
  G4double GetSurfaceArea() override;
 
  G4GeometryType GetEntityType() const override;
 
  G4VisExtent GetExtent() const override;
 
  void BoundingLimits(G4ThreeVector& pMin, G4ThreeVector& pMax) const override;
 
  G4bool CalculateExtent(const EAxis pAxis, const G4VoxelLimits& pVoxelLimit,
                         const G4AffineTransform& pTransform, G4double& pMin,
                         G4double& pMax) const override;
 
  void DescribeYourselfTo(G4VGraphicsScene& scene) const override;
 
  G4Polyhedron* CreatePolyhedron() const override;
 
  std::ostream& StreamInfo(std::ostream& os) const override;
 
  // ── G4OCCTSolid-specific interface ────────────────────────────────────────
 
  const TopoDS_Shape& GetOCCTShape() const { return fShape; }
 
  void SetOCCTShape(const TopoDS_Shape& shape) {
    if (shape.IsNull()) {
      throw std::invalid_argument("G4OCCTSolid::SetOCCTShape: shape must not be null");
    }
    fShape = shape;
    ComputeBounds();
    {
      std::unique_lock<std::mutex> lock(fPolyhedronMutex);
      fCachedPolyhedron.reset();
    }
    {
      std::unique_lock<std::mutex> lock(fVolumeAreaMutex);
      fCachedVolume.reset();
      fCachedSurfaceArea.reset();
    }
    {
      std::unique_lock<std::mutex> lock(fSurfaceCacheMutex);
      fSurfaceCache.reset();
    }
    fShapeGeneration.fetch_add(1, std::memory_order_release);
  }
 
private:
  struct FaceBounds {
    TopoDS_Face face;
    Bnd_Box box;
    BRepAdaptor_Surface adaptor;
    std::optional<gp_Pln> plane; 
    std::vector<gp_Pnt2d> uvPolygon;
    std::optional<G4ThreeVector> outwardNormal;
  };
 
  struct ClosestFaceMatch {
    TopoDS_Face face;
    G4double distance{kInfinity};
    std::size_t faceIndex{0};
    std::optional<std::pair<Standard_Real, Standard_Real>> uv;
  };
 
  struct AxisAlignedBounds {
    G4ThreeVector min;
    G4ThreeVector max;
  };
 
  struct ClassifierCache {
    std::uint64_t generation{std::numeric_limits<std::uint64_t>::max()};
    std::optional<BRepClass3d_SolidClassifier> classifier;
  };
 
  struct IntersectorCache {
    std::uint64_t generation{std::numeric_limits<std::uint64_t>::max()};
    std::vector<std::unique_ptr<IntCurvesFace_Intersector>> faceIntersectors;
    std::vector<Bnd_Box> expandedBoxes; 
  };
 
  struct InscribedSphere {
    G4ThreeVector centre;
    G4double radius;
  };
 
  struct SphereCacheData {
    std::vector<InscribedSphere> spheres;
    std::uint64_t generation{std::numeric_limits<std::uint64_t>::max()};
  };
 
  static constexpr std::size_t kMaxInscribedSpheres = 64;
 
  struct SurfaceTriangle {
    G4ThreeVector p1, p2, p3;
    std::uint32_t
        faceIndex; 
  };
 
  struct SurfaceSamplingCache {
    std::vector<TopoDS_Face> faces; 
    std::vector<SurfaceTriangle> triangles;
    std::vector<G4double> cumulativeAreas;
    G4double totalArea{0.0};
  };
 
  TopoDS_Shape fShape;
 
  AxisAlignedBounds fCachedBounds;
 
  std::vector<FaceBounds> fFaceBoundsCache;
 
  std::unordered_map<const TopoDS_TShape*, std::vector<std::size_t>> fFaceAdaptorIndex;
 
  std::vector<InscribedSphere> fInitialSpheres;
 
  Handle(BRepExtrema_TriangleSet) fTriangleSet;
 
  G4double fBVHDeflection{0.0};
 
  std::vector<G4double> fFaceDeflections;
 
  bool fAllFacesPlanar{false};
 
  std::atomic<std::uint64_t> fShapeGeneration{0};
 
  mutable G4Cache<ClassifierCache> fClassifierCache;
 
  mutable G4Cache<IntersectorCache> fIntersectorCache;
 
  mutable G4Cache<SphereCacheData> fSphereCache;
 
  mutable std::unique_ptr<G4Polyhedron> fCachedPolyhedron;
 
  mutable std::uint64_t fPolyhedronGeneration{std::numeric_limits<std::uint64_t>::max()};
 
  mutable bool fPolyhedronBuilding{false};
 
  mutable std::mutex fPolyhedronMutex;
 
  mutable std::condition_variable fPolyhedronCV;
 
  mutable std::optional<G4double> fCachedVolume;
 
  mutable std::optional<G4double> fCachedSurfaceArea;
 
  mutable std::mutex fVolumeAreaMutex;
 
  mutable std::optional<SurfaceSamplingCache> fSurfaceCache;
 
  mutable std::uint64_t fSurfaceCacheGeneration{std::numeric_limits<std::uint64_t>::max()};
 
  mutable std::mutex fSurfaceCacheMutex;
 
  BRepClass3d_SolidClassifier& GetOrCreateClassifier() const;
 
  IntersectorCache& GetOrCreateIntersector() const;
 
  SphereCacheData& GetOrInitSphereCache() const;
 
  void TryInsertSphere(const G4ThreeVector& centre, G4double d) const;
 
  void ComputeBounds();
 
  void ComputeInitialSpheres();
 
  G4double AABBLowerBound(const G4ThreeVector& p) const;
 
  G4double BVHLowerBoundDistance(const G4ThreeVector& p) const;
 
  G4double PlanarFaceLowerBoundDistance(const G4ThreeVector& p) const;
 
  static std::optional<ClosestFaceMatch>
  TryFindClosestFace(const std::vector<FaceBounds>& faceBoundsCache, const G4ThreeVector& point,
                     G4double maxDistance = kInfinity);
 
  const SurfaceSamplingCache& GetOrBuildSurfaceCache() const;
};
 
#endif // G4OCCT_G4OCCTSolid_hh