From 11ecb3a0d9bfd6ed0942305e868dfc4cc286a7f5 Mon Sep 17 00:00:00 2001
From: Evgueni Tcherniaev <evgueni.tcherniaev@cern.ch>
Date: Tue, 27 Aug 2024 19:21:32 +0200
Subject: [PATCH] Extended functionality related to volume calculation

---
 .../GeoModelKernel/GeoModelKernel/GeoBox.h    | 10 +++-
 .../GeoModelKernel/GeoModelKernel/GeoCons.h   | 10 ++--
 .../GeoModelKernel/GeoEllipticalTube.h        |  8 ++-
 .../GeoModelKernel/GeoGenericTrap.h           |  7 ++-
 .../GeoModelKernel/GeoModelKernel/GeoPara.h   |  7 ++-
 .../GeoModelKernel/GeoModelKernel/GeoPcon.h   |  9 ++-
 .../GeoModelKernel/GeoModelKernel/GeoPgon.h   |  8 ++-
 .../GeoModelKernel/GeoModelKernel/GeoShape.h  | 28 ++++++++--
 .../GeoModelKernel/GeoShapeIntersection.h     | 15 +++--
 .../GeoModelKernel/GeoShapeShift.h            | 16 ++++--
 .../GeoModelKernel/GeoShapeSubtraction.h      | 15 +++--
 .../GeoModelKernel/GeoShapeUnion.h            | 21 ++++---
 .../GeoModelKernel/GeoSimplePolygonBrep.h     |  7 ++-
 .../GeoModelKernel/GeoTessellatedSolid.h      |  7 ++-
 .../GeoModelKernel/GeoModelKernel/GeoTorus.h  |  9 ++-
 .../GeoModelKernel/GeoModelKernel/GeoTrap.h   |  8 ++-
 .../GeoModelKernel/GeoModelKernel/GeoTrd.h    | 11 +++-
 .../GeoModelKernel/GeoModelKernel/GeoTube.h   | 11 +++-
 .../GeoModelKernel/GeoModelKernel/GeoTubs.h   |  7 ++-
 .../GeoModelKernel/GeoTwistedTrap.h           |  8 ++-
 .../GeoModelKernel/GeoUnidentifiedShape.h     |  5 ++
 GeoModelCore/GeoModelKernel/src/GeoBox.cxx    |  7 +--
 GeoModelCore/GeoModelKernel/src/GeoCons.cxx   |  8 +--
 .../GeoModelKernel/src/GeoEllipticalTube.cxx  | 11 ++--
 .../GeoModelKernel/src/GeoGenericTrap.cxx     |  6 +-
 GeoModelCore/GeoModelKernel/src/GeoPara.cxx   |  4 +-
 GeoModelCore/GeoModelKernel/src/GeoPcon.cxx   |  6 +-
 GeoModelCore/GeoModelKernel/src/GeoPgon.cxx   |  5 +-
 GeoModelCore/GeoModelKernel/src/GeoShape.cxx  |  8 +--
 .../src/GeoShapeIntersection.cxx              | 31 ++++++-----
 .../GeoModelKernel/src/GeoShapeShift.cxx      | 19 +++----
 .../src/GeoShapeSubtraction.cxx               | 55 +++++++++++++------
 .../GeoModelKernel/src/GeoShapeUnion.cxx      | 51 ++++++++++++-----
 .../src/GeoSimplePolygonBrep.cxx              |  8 +--
 .../src/GeoTessellatedSolid.cxx               |  2 +-
 GeoModelCore/GeoModelKernel/src/GeoTorus.cxx  |  6 +-
 GeoModelCore/GeoModelKernel/src/GeoTrap.cxx   |  7 +--
 GeoModelCore/GeoModelKernel/src/GeoTrd.cxx    |  3 +-
 GeoModelCore/GeoModelKernel/src/GeoTube.cxx   |  9 +--
 GeoModelCore/GeoModelKernel/src/GeoTubs.cxx   |  7 +--
 .../GeoModelKernel/src/GeoTwistedTrap.cxx     |  2 +-
 41 files changed, 306 insertions(+), 176 deletions(-)

diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoBox.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoBox.h
index 3542adc22..d823885ac 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoBox.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoBox.h
@@ -14,7 +14,7 @@ class GeoBox : public GeoShape
   GeoBox (double XHalfLength, double YHalfLength, double ZHalfLength);
 
   //    Returns the volume of the shape, for mass inventory
-  virtual double volume () const;
+  virtual double volume (int npoints = 0) const;
 
   //    Returns the bonding box of the shape
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -33,7 +33,12 @@ class GeoBox : public GeoShape
      return getClassTypeID();
   }
 
-  //     Executes a GeoShapeAction
+  //    Returns true as BOX is a polyhedron.
+  virtual bool isPolyhedron () const {
+    return true;
+  }
+
+  //    Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const;
 
   //    For type identification
@@ -61,7 +66,6 @@ class GeoBox : public GeoShape
      return m_zHalfLength;
   }
 
-
  protected:
   virtual ~GeoBox() = default;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoCons.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoCons.h
index b7fbf473c..722131221 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoCons.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoCons.h
@@ -12,7 +12,7 @@ class GeoCons : public GeoShape {
   GeoCons (double RMin1, double RMin2, double RMax1, double RMax2, double DZ, double SPhi, double DPhi);
 
   //    Returns the volume of the shape, for mass inventory.
-  virtual double volume () const;
+  virtual double volume (int npoints = 0) const;
 
   //    Returns the bonding box of the shape.
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -31,6 +31,11 @@ class GeoCons : public GeoShape {
      return getClassTypeID();
   }
 
+  //    Returns false as CONS is not a polyhedron.
+  virtual bool isPolyhedron () const {
+    return false;
+  }
+
   //    Executes a GeoShapeAction.
   virtual void exec (GeoShapeAction *action) const;
 
@@ -63,8 +68,6 @@ protected:
   virtual ~GeoCons() = default;
 
  private:
-
-
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
@@ -77,5 +80,4 @@ protected:
   double m_dPhi{0.};
 };
 
-
 #endif
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoEllipticalTube.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoEllipticalTube.h
index a36c03120..79dfdb04c 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoEllipticalTube.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoEllipticalTube.h
@@ -20,7 +20,7 @@ class GeoEllipticalTube : public GeoShape {
   GeoEllipticalTube(double XHalfLength, double YHalfLength, double ZHalfLength);
 
   //    Returns the volume of the shape, for mass inventory
-  virtual double volume () const;
+  virtual double volume (int npoints = 0) const;
 
   //    Returns the bonding box of the shape
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -39,6 +39,11 @@ class GeoEllipticalTube : public GeoShape {
     return getClassTypeID();
   }
 
+  //    Returns false as ELLIPTICAL TUBE is not a polyhedron.
+  virtual bool isPolyhedron () const {
+    return false;
+  }
+
   //    Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const;
 
@@ -71,7 +76,6 @@ class GeoEllipticalTube : public GeoShape {
   virtual ~GeoEllipticalTube() = default;
 
  private:
-
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoGenericTrap.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoGenericTrap.h
index edbc1c077..4cf59192d 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoGenericTrap.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoGenericTrap.h
@@ -17,7 +17,7 @@ class GeoGenericTrap : public GeoShape {
   GeoGenericTrap(double ZHalfLength, const GeoGenericTrapVertices& Vertices);
 
   //    Returns the volume of the shape, for mass inventory.
-  virtual double volume() const;
+  virtual double volume(int npoints = 0) const;
 
   //    Returns the bonding box of the shape.
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -36,6 +36,11 @@ class GeoGenericTrap : public GeoShape {
      return getClassTypeID();
   }
 
+  //    Returns false as GENERIC TRAP is not a polyhedron.
+  virtual bool isPolyhedron () const {
+    return false;
+  }
+
   //    For type identification.
   static const std::string& getClassType() {
     return s_classType;
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPara.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPara.h
index 26ad8555e..c98a41f15 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPara.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPara.h
@@ -14,7 +14,7 @@ class GeoPara : public GeoShape
   GeoPara (double XHalfLength, double YHalfLength, double ZHalfLength, double Alpha, double Theta, double Phi);
 
   //    Returns the volume of the shape, for mass inventory.
-  virtual double volume () const;
+  virtual double volume (int npoints = 0) const;
 
   //    Returns the bonding box of the shape.
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -33,6 +33,11 @@ class GeoPara : public GeoShape
     return getClassTypeID();
   }
 
+  //    Returns true as PARA is a polyhedron.
+  virtual bool isPolyhedron () const {
+    return true;
+  }
+
   //    Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPcon.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPcon.h
index 59b88f989..080fd1ecf 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPcon.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPcon.h
@@ -25,7 +25,7 @@ class GeoPcon : public GeoShape {
   GeoPcon (double SPhi, double DPhi);
 
   //    Returns the volume of the shape, for mass inventory
-  virtual double volume () const;
+  virtual double volume (int npoints = 0) const;
 
   //    Returns the bonding box of the shape
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -44,6 +44,11 @@ class GeoPcon : public GeoShape {
      return getClassTypeID();
   }
 
+  //    Returns false as PCON is not a polyhedron.
+  virtual bool isPolyhedron () const {
+    return false;
+  }
+
   //    Add another plane to the polycone  A minimum of two
   //    planes are required to create a valid polycone.
   void addPlane (double ZPlane, double RMinPlane, double RMaxPlane);
@@ -102,7 +107,6 @@ class GeoPcon : public GeoShape {
   virtual ~GeoPcon() = default;
 
  private:
-
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
@@ -119,5 +123,4 @@ class GeoPcon : public GeoShape {
   std::vector<double> m_rMaxPlane{};
 };
 
-
 #endif
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPgon.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPgon.h
index 7dea7fad2..4ca4e1d9e 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPgon.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPgon.h
@@ -17,7 +17,7 @@ class GeoPgon : public GeoShape
   GeoPgon (double SPhi, double DPhi, unsigned int NSides);
 
   //    Returns the volume of the shape, for mass inventory
-  virtual double volume () const;
+  virtual double volume (int npoints = 0) const;
 
   //    Returns the bonding box of the shape
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -36,6 +36,11 @@ class GeoPgon : public GeoShape
      return getClassTypeID();
   }
 
+  //    Returns true as PGON is a polyhedron.
+  virtual bool isPolyhedron () const {
+    return true;
+  }
+
   //    Add another plane to the polygon. A minimum of two
   //    planes are required to create a valid polygon.
   void addPlane (double ZPlane, double RMinPlane, double RMaxPlane);
@@ -99,7 +104,6 @@ class GeoPgon : public GeoShape
   virtual ~GeoPgon() = default;
 
  private:
-  
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShape.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShape.h
index 493ef7046..452139ca1 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShape.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShape.h
@@ -27,6 +27,7 @@
 #include "GeoModelKernel/RCBase.h"
 #include <GeoModelKernel/GeoDefinitions.h>
 #include <string>
+#include <atomic>
 
 using ShapeType = unsigned int;
 class GeoShapeIntersection;
@@ -42,7 +43,7 @@ class GeoShape : public RCBase
   GeoShape () = default;
 
   //    Returns the volume of the shape, for mass inventory.
-  virtual double volume () const;
+  virtual double volume (int npoints = 1000000) const;
 
   //    Returns the bonding box of the shape.
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -69,19 +70,38 @@ class GeoShape : public RCBase
   //    Returns the shape type, as an coded integer.
   virtual ShapeType typeID () const = 0;
 
+  //    Returns true if the shape is a polyhedron.
+  virtual bool isPolyhedron () const = 0;
+
+  //    Returns number of constituents
+  virtual unsigned int getNoConstituents () const {
+    return 1;
+  }
+
   //    Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const = 0;
 
+  //    Returns volume value
+  double getVolumeValue () const {
+    return m_shape_volume.load();
+  }
+
+  //    Sets volume value
+  void setVolumeValue (double value) const {
+    m_shape_volume = value;
+  }
+
  protected:
-  virtual ~GeoShape() = default;
+  virtual ~GeoShape () = default;
 
   //    Returns the bounding box of the specified disk. This method is used
   //    for calculation of the extend of a tube, cone, polycone, torus.
   static void diskExtent(double rmin, double rmax, double sphi, double dphi,
                          double& xmin, double& ymin, double& xmax, double& ymax);
 
-
-
+ private:
+  //    Cached volume
+  mutable std::atomic<double> m_shape_volume{-1.};
 };
 
 #endif
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeIntersection.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeIntersection.h
index 7c1027339..bb8f4dd2a 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeIntersection.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeIntersection.h
@@ -17,7 +17,7 @@ class GeoShapeIntersection : public GeoShape
   GeoShapeIntersection (const GeoShape* A, const GeoShape* B);
 
   //    Returns the volume of the shape, for mass inventory
-  virtual double volume () const;
+  virtual double volume (int npoints = 1000000) const;
 
   //    Returns the bonding box of the shape
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -36,6 +36,16 @@ class GeoShapeIntersection : public GeoShape
      return getClassTypeID();
   }
 
+  //    Returns true if the resulting shape is a polyhedron, false otherwise.
+  virtual bool isPolyhedron () const {
+    return (m_opA->isPolyhedron() && m_opB->isPolyhedron());
+  }
+
+  //    Returns total number of constituents
+  virtual unsigned int getNoConstituents () const {
+    return (m_opA->getNoConstituents() + m_opB->getNoConstituents());
+  }
+
   //    Returns the first operand being ANDed
   const GeoShape* getOpA () const {
       return m_opA;
@@ -68,9 +78,6 @@ class GeoShapeIntersection : public GeoShape
   //    The second shape operand in the AND operation.
   GeoIntrusivePtr<const GeoShape> m_opB{};
 
-  //    Cached volume
-  mutable std::atomic<double> fVolume{-1};
-
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeShift.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeShift.h
index ab7fa80b8..1ce417748 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeShift.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeShift.h
@@ -17,8 +17,8 @@ class GeoShapeShift : public GeoShape {
   GeoShapeShift (const GeoShape* A, const GeoTrf::Transform3D &X);
 
   //    Returns the volume of the shape, for mass inventory
-  virtual double volume () const {
-     return m_op->volume();
+  virtual double volume (int npoints = 1000000) const {
+     return m_op->volume(npoints);
   }
 
   //    Returns the bonding box of the shape
@@ -38,6 +38,16 @@ class GeoShapeShift : public GeoShape {
     return getClassTypeID();
   }
 
+  //    Returns true if OR is a polyhedron, false otherwise
+  virtual bool isPolyhedron () const {
+    return m_op->isPolyhedron();
+  }
+
+  //    Returns number of constituents
+  virtual unsigned int getNoConstituents () const {
+    return m_op->getNoConstituents();
+  }
+
   // Returns the first operand being ORed
   const GeoShape* getOp() const {
       return m_op;
@@ -51,13 +61,11 @@ class GeoShapeShift : public GeoShape {
   //    Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const;
 
- 
 
   static ShapeType getClassTypeID () {
      return s_classTypeID;
   }
 
-
   //    For type identification.
   static const std::string& getClassType () {
      return s_classType;
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeSubtraction.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeSubtraction.h
index 6036e0921..ca9cb5816 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeSubtraction.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeSubtraction.h
@@ -16,7 +16,7 @@ class GeoShapeSubtraction : public GeoShape {
   GeoShapeSubtraction (const GeoShape* A, const GeoShape* B);
 
   //    Returns the volume of the shape, for mass inventory
-  virtual double volume () const;
+  virtual double volume (int npoints = 1000000) const;
 
   //    Returns the bonding box of the shape
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -35,6 +35,16 @@ class GeoShapeSubtraction : public GeoShape {
      return  getClassTypeID();
   }
 
+  //    Returns true if the resulting shape is a polyhedron, false otherwise.
+  virtual bool isPolyhedron () const {
+    return (m_opA->isPolyhedron() && m_opB->isPolyhedron());
+  }
+
+  //    Returns total number of constituents
+  virtual unsigned int getNoConstituents () const {
+    return (m_opA->getNoConstituents() + m_opB->getNoConstituents());
+  }
+
   //    Returns the first operand in the subtraction
   const GeoShape* getOpA () const{
      return m_opA;
@@ -69,9 +79,6 @@ class GeoShapeSubtraction : public GeoShape {
   //    The second shape operand in the Subtraction operation
   GeoIntrusivePtr<const GeoShape> m_opB{nullptr};
 
-  //    Cached volume
-  mutable std::atomic<double> fVolume{-1.};
-
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeUnion.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeUnion.h
index 7849753e5..47b66b8ce 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeUnion.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeUnion.h
@@ -16,7 +16,7 @@ class GeoShapeUnion : public GeoShape
   GeoShapeUnion (const GeoShape* A, const GeoShape* B);
 
   //    Returns the volume of the shape, for mass inventory
-  virtual double volume () const;
+  virtual double volume (int npoints = 1000000) const;
 
   //    Returns the bonding box of the shape
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -24,8 +24,6 @@ class GeoShapeUnion : public GeoShape
 
   //    Returns true if the shape contains the point, false otherwise
   virtual bool contains (double x, double y, double z) const;
-
- 
  
   //    Returns the OR shape type, as a string.
   virtual const std::string & type() const{
@@ -37,6 +35,16 @@ class GeoShapeUnion : public GeoShape
      return getClassTypeID();
   }
 
+  //    Returns true if the resulting shape is a polyhedron, false otherwise.
+  virtual bool isPolyhedron () const {
+    return (m_opA->isPolyhedron() && m_opB->isPolyhedron());
+  }
+
+  //    Returns total number of constituents
+  virtual unsigned int getNoConstituents () const {
+    return (m_opA->getNoConstituents() + m_opB->getNoConstituents());
+  }
+
   //    Returns the first operand being ORed
   const GeoShape* getOpA() const {
       return m_opA;
@@ -64,21 +72,16 @@ class GeoShapeUnion : public GeoShape
   virtual ~GeoShapeUnion() = default;
 
  private:
-
-
   //    The first shape operand in the OR operation
   GeoIntrusivePtr<const GeoShape> m_opA{};
 
   //    The second shape operand in the OR operation
   GeoIntrusivePtr<const GeoShape> m_opB{};
 
-  //    Cached volume
-  mutable std::atomic<double> fVolume{-1.};
-
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
-    //    For I/O only!
+  //    For I/O only!
   GeoShapeUnion() = default;
   friend Persistifier;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoSimplePolygonBrep.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoSimplePolygonBrep.h
index 6a1488091..8e55f021e 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoSimplePolygonBrep.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoSimplePolygonBrep.h
@@ -27,7 +27,7 @@ class GeoSimplePolygonBrep : public GeoShape {
   GeoSimplePolygonBrep(double dz);
 
   //    Returns the volume of the shape, for mass inventory
-  virtual double volume() const;
+  virtual double volume (int npoints = 0) const;
 
   //    Returns the bonding box of the shape
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -46,6 +46,11 @@ class GeoSimplePolygonBrep : public GeoShape {
     return getClassTypeID();
   }
 
+  //    Returns true as BREP is a polyhedron.
+  virtual bool isPolyhedron () const {
+    return true;
+  }
+
   //    For type identification
   static const std::string& getClassType() {
     return s_classType;
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTessellatedSolid.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTessellatedSolid.h
index 72850483c..44980db31 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTessellatedSolid.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTessellatedSolid.h
@@ -17,7 +17,7 @@ class GeoTessellatedSolid : public GeoShape {
   GeoTessellatedSolid();
 
   //    Returns the volume of the shape, for mass inventory
-  virtual double volume() const;
+  virtual double volume(int npoints = 0) const;
 
   //    Returns the bonding box of the shape
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -36,6 +36,11 @@ class GeoTessellatedSolid : public GeoShape {
      return getClassTypeID();
   }
 
+  //    Returns true as TESSELLATED SOLID is a polyhedron.
+  virtual bool isPolyhedron () const {
+    return true;
+  }
+
   //    For type identification
   static const std::string& getClassType();
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTorus.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTorus.h
index cbbe6fd2b..2dba24c3a 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTorus.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTorus.h
@@ -13,7 +13,7 @@ class GeoTorus : public GeoShape {
   GeoTorus (double Rmin, double Rmax, double Rtor, double SPhi, double DPhi);
 
   //    Returns the volume of the shape, for mass inventory.
-  virtual double volume () const;
+  virtual double volume (int npoints = 0) const;
 
   //    Returns the bonding box of the shape.
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -32,6 +32,11 @@ class GeoTorus : public GeoShape {
      return getClassTypeID();
   }
 
+  //    Returns false as TORUS is not a polyhedron.
+  virtual bool isPolyhedron () const {
+    return false;
+  }
+
   //    Executes a GeoShapeAction.
   virtual void exec (GeoShapeAction *action) const;
 
@@ -74,8 +79,6 @@ class GeoTorus : public GeoShape {
   virtual ~GeoTorus() = default;
 
  private:
-
-
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrap.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrap.h
index 3726befa6..22b4cb60c 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrap.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrap.h
@@ -19,7 +19,7 @@ class GeoTrap : public GeoShape
   GeoTrap (double ZHalfLength, double Theta, double Phi, double Dydzn, double Dxdyndzn, double Dxdypdzn, double Angleydzn, double Dydzp, double Dxdyndzp, double Dxdypdzp, double Angleydzp);
 
   //    Returns the volume of the shape, for mass inventory.
-  virtual double volume () const;
+  virtual double volume (int npoints = 0) const;
 
   //    Returns the bonding box of the shape.
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -38,6 +38,11 @@ class GeoTrap : public GeoShape
      return getClassTypeID();
   }
 
+  //    Returns true as TRAP is a polyhedron.
+  virtual bool isPolyhedron () const {
+    return true;
+  }
+
   //    Executes a GeoShapeAction.
   virtual void exec (GeoShapeAction *action) const;
 
@@ -86,7 +91,6 @@ class GeoTrap : public GeoShape
 
  private:
 
-
   static const std::string s_classType;
   static const ShapeType s_classTypeID;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrd.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrd.h
index 48b6759e5..1912351c6 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrd.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrd.h
@@ -13,7 +13,7 @@ class GeoTrd : public GeoShape
   GeoTrd (double XHalfLength1, double XHalfLength2, double YHalfLength1, double YHalfLength2, double ZHalfLength);
 
   //    Returns the volume of the shape, for mass inventory.
-  virtual double volume () const;
+  virtual double volume (int npoints = 0) const;
 
   //    Returns the bonding box of the shape.
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -32,6 +32,11 @@ class GeoTrd : public GeoShape
      return getClassTypeID();
   }
 
+  //    Returns true as TRD is a polyhedron.
+  virtual bool isPolyhedron () const {
+    return true;
+  }
+
   //    Executes a GeoShapeAction.
   virtual void exec (GeoShapeAction *action) const;
 
@@ -51,10 +56,10 @@ class GeoTrd : public GeoShape
   //    Half length in the x-direction at +dz.
   double  getXHalfLength2 () const { return m_xHalfLength2; }
 
-  //    Half-length in the y direction at +dz.
+  //    Half-length in the y direction at -dz.
   double  getYHalfLength1 () const { return m_yHalfLength1; } 
 
-  //    Half-length in the y direction at -dz.
+  //    Half-length in the y direction at +dz.
   double  getYHalfLength2 () const { return m_yHalfLength2; } 
 
   //    Half-length in the z direction.
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTube.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTube.h
index f9d7e2370..46630438d 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTube.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTube.h
@@ -13,7 +13,7 @@ class GeoTube : public GeoShape
   GeoTube (double RMin, double RMax, double ZHalfLength);
 
   //    Returns the volume of the shape, for mass inventory.
-  virtual double volume () const;
+  virtual double volume (int npoints = 0) const;
 
   //    Returns the bonding box of the shape.
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -31,6 +31,12 @@ class GeoTube : public GeoShape
   virtual ShapeType typeID () const{
      return getClassTypeID();
   }
+
+  //    Returns false as TUBE is not a polyhedron.
+  virtual bool isPolyhedron () const {
+    return false;
+  }
+
   //    Executes a GeoShapeAction.
   virtual void exec (GeoShapeAction *action) const;
 
@@ -60,8 +66,7 @@ class GeoTube : public GeoShape
   }
 
  protected:
-  //## Destructor (generated)
-  virtual ~GeoTube();
+  virtual ~GeoTube() = default;
 
  private:
   GeoTube(const GeoTube &right);
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTubs.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTubs.h
index 79d3dbf40..5e41b2238 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTubs.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTubs.h
@@ -13,7 +13,7 @@ class GeoTubs : public GeoShape
   GeoTubs (double RMin, double RMax, double ZHalfLength, double SPhi, double DPhi);
 
   //    Returns the volume of the shape, for mass inventory.
-  virtual double volume () const;
+  virtual double volume (int npoints = 0) const;
 
   //    Returns the bonding box of the shape.
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -32,6 +32,11 @@ class GeoTubs : public GeoShape
      return getClassTypeID();
   }
 
+  //    Returns false as TUBS is not a polyhedron.
+  virtual bool isPolyhedron () const {
+    return false;
+  }
+
   //    Executes a GeoShapeAction.
   virtual void exec (GeoShapeAction *action) const;
 
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTwistedTrap.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTwistedTrap.h
index e54dd5293..6142d633e 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTwistedTrap.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTwistedTrap.h
@@ -64,7 +64,7 @@ class GeoTwistedTrap : public GeoShape
   );
 
   //    Returns the volume of the shape, for mass inventory
-  virtual double volume () const;
+  virtual double volume (int npoints = 0) const;
 
   //    Returns the bonding box of the shape.
   virtual void extent (double& xmin, double& ymin, double& zmin,
@@ -83,6 +83,11 @@ class GeoTwistedTrap : public GeoShape
      return getClassTypeID();
   }
 
+  //    Returns false as TWISTED TRAP is not a polyhedron.
+  virtual bool isPolyhedron () const {
+    return false;
+  }
+
   //    Executes a GeoShapeAction
   virtual void exec (GeoShapeAction *action) const;
 
@@ -96,7 +101,6 @@ class GeoTwistedTrap : public GeoShape
      return s_classTypeID;
   }
 
-
   inline double getY1HalfLength() const { return m_dy1 ; }
   inline double getX1HalfLength() const { return m_dx1 ; }
   inline double getX2HalfLength() const { return m_dx2 ; }
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoUnidentifiedShape.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoUnidentifiedShape.h
index 8f021aee8..34d4efa1e 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoUnidentifiedShape.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoUnidentifiedShape.h
@@ -54,6 +54,11 @@ class GeoUnidentifiedShape: public GeoShape {
   // Returns the shape type, as an coded integer.
   virtual ShapeType typeID () const;
 
+  // Returns false as the shape is not a polyhedron.
+  virtual bool isPolyhedron () const {
+    return false;
+  }
+
   // For type identification.
   static const std::string& getClassType ();
 
diff --git a/GeoModelCore/GeoModelKernel/src/GeoBox.cxx b/GeoModelCore/GeoModelKernel/src/GeoBox.cxx
index 209d007c6..327b72cee 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoBox.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoBox.cxx
@@ -11,10 +11,10 @@ const ShapeType GeoBox::s_classTypeID = 0x10; // 16
 GeoBox::GeoBox (double XHalfLength, double YHalfLength, double ZHalfLength):
   m_xHalfLength {XHalfLength},
   m_yHalfLength {YHalfLength},
-  m_zHalfLength {ZHalfLength} {}
+  m_zHalfLength {ZHalfLength}
+{}
 
-
-double GeoBox::volume () const {
+double GeoBox::volume (int) const {
   return 8.0 * m_xHalfLength * m_yHalfLength * m_zHalfLength;
 }
 
@@ -35,7 +35,6 @@ bool GeoBox::contains (double x, double y, double z) const {
   return (std::max(std::max(distx, disty), distz) <= 0.0);
 }
 
-
 void GeoBox::exec (GeoShapeAction *action) const {
   action->handleBox(this);
 }
diff --git a/GeoModelCore/GeoModelKernel/src/GeoCons.cxx b/GeoModelCore/GeoModelKernel/src/GeoCons.cxx
index 2fcad24bf..4c0809462 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoCons.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoCons.cxx
@@ -17,11 +17,9 @@ GeoCons::GeoCons (double RMin1, double RMin2, double RMax1, double RMax2, double
   , m_dZ (DZ)
   , m_sPhi (SPhi)
   , m_dPhi (DPhi)
-{
-}
-
+{}
 
-double GeoCons::volume () const
+double GeoCons::volume (int) const
 {
   return (m_dPhi * (1./3.)) * m_dZ *
     (m_rMax1 * m_rMax1 + m_rMax2 * m_rMax2 + m_rMax1 * m_rMax2 -
@@ -59,8 +57,6 @@ bool GeoCons::contains (double x, double y, double z) const
   return (m_dPhi <= M_PI) ? (ds <= 0 && de <= 0) : (ds <= 0 || de <= 0);
 }
 
-
-
 void GeoCons::exec (GeoShapeAction *action) const {
   action->handleCons(this);
 }
diff --git a/GeoModelCore/GeoModelKernel/src/GeoEllipticalTube.cxx b/GeoModelCore/GeoModelKernel/src/GeoEllipticalTube.cxx
index e356b15e1..f1e9e487d 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoEllipticalTube.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoEllipticalTube.cxx
@@ -11,12 +11,12 @@ const std::string GeoEllipticalTube::s_classType = "EllipticalTube";
 const ShapeType GeoEllipticalTube::s_classTypeID = 0x22;
 
 GeoEllipticalTube::GeoEllipticalTube(double XHalfLength, double YHalfLength, double ZHalfLength): 
-    m_xHalfLength{XHalfLength},
-   m_yHalfLength{YHalfLength},
-   m_zHalfLength{ZHalfLength} {}
+  m_xHalfLength{XHalfLength},
+  m_yHalfLength{YHalfLength},
+  m_zHalfLength{ZHalfLength}
+{}
 
-
-double GeoEllipticalTube::volume () const
+double GeoEllipticalTube::volume (int) const
 {
 #ifndef M_PI
   constexpr double M_PI = 3.14159265358979323846;
@@ -42,7 +42,6 @@ bool GeoEllipticalTube::contains (double x, double y, double z) const
           (y * y) / (m_yHalfLength * m_yHalfLength) <= 1.0);
 }
 
-
 void GeoEllipticalTube::exec (GeoShapeAction *action) const {
   action->handleEllipticalTube(this);
 }
diff --git a/GeoModelCore/GeoModelKernel/src/GeoGenericTrap.cxx b/GeoModelCore/GeoModelKernel/src/GeoGenericTrap.cxx
index 31e1b31b8..fe87f7fc0 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoGenericTrap.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoGenericTrap.cxx
@@ -10,10 +10,10 @@ const ShapeType GeoGenericTrap::s_classTypeID = 0x23;
 
 GeoGenericTrap::GeoGenericTrap(double ZHalfLength, const GeoGenericTrapVertices& Vertices)
   : m_zHalfLength(ZHalfLength)
-  , m_vertices(Vertices) {}
+  , m_vertices(Vertices)
+{}
 
-
-double GeoGenericTrap::volume() const
+double GeoGenericTrap::volume(int) const
 {
   // diagonals
   GeoTwoVector a = m_vertices[3] - m_vertices[1];
diff --git a/GeoModelCore/GeoModelKernel/src/GeoPara.cxx b/GeoModelCore/GeoModelKernel/src/GeoPara.cxx
index c1634dee1..cd195a954 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoPara.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoPara.cxx
@@ -19,9 +19,7 @@ m_alpha (Alpha),
 m_phi (Phi)
 {}
 
-
-
-double GeoPara::volume () const
+double GeoPara::volume (int) const
 {
   return 8.0 * m_xHalfLength * m_yHalfLength * m_zHalfLength;
 }
diff --git a/GeoModelCore/GeoModelKernel/src/GeoPcon.cxx b/GeoModelCore/GeoModelKernel/src/GeoPcon.cxx
index 858cfdf12..23c617cb9 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoPcon.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoPcon.cxx
@@ -13,11 +13,9 @@ const ShapeType GeoPcon::s_classTypeID = 0x13;
 GeoPcon::GeoPcon (double SPhi, double DPhi)
   : m_sPhi (SPhi)
   , m_dPhi (DPhi)
-{
-}
-
+{}
 
-double GeoPcon::volume () const
+double GeoPcon::volume (int) const
 {
   if (!isValid ())
     throw std::runtime_error ("Volume requested for incomplete polycone");
diff --git a/GeoModelCore/GeoModelKernel/src/GeoPgon.cxx b/GeoModelCore/GeoModelKernel/src/GeoPgon.cxx
index dd8a2d9e6..2ad286208 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoPgon.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoPgon.cxx
@@ -14,10 +14,9 @@ GeoPgon::GeoPgon (double SPhi, double DPhi, unsigned int NSides)
   : m_sPhi (SPhi)
   , m_dPhi (DPhi)
   , m_nSides (NSides)
-{
-}
+{}
 
-double GeoPgon::volume () const
+double GeoPgon::volume (int) const
 {
 #ifndef M_PI
   constexpr double M_PI = 3.14159265358979323846;
diff --git a/GeoModelCore/GeoModelKernel/src/GeoShape.cxx b/GeoModelCore/GeoModelKernel/src/GeoShape.cxx
index 9049e957b..78aa1cb78 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoShape.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoShape.cxx
@@ -7,15 +7,15 @@
 #include "GeoModelKernel/GeoShapeUnion.h"
 #include "GeoModelKernel/GeoShapeSubtraction.h"
 #include "GeoModelKernel/GeoShapeShift.h"
+#include <algorithm>
 #include <cmath>
 #include <cstdint>
 
-
-double GeoShape::volume () const
+double GeoShape::volume (int npoints) const
 {
-  constexpr int npoints = 1000000;     // number of random points
   constexpr double expansion = 0.001;  // bounding box expansion
   constexpr double f = 1./4294967296.; // 2^-32 - int to double conversion
+  int np = std::max(npoints, 1000);    // number of points is at least 1000
 
   // set up bonding box
   double xmin = 0, ymin = 0, zmin = 0, xmax = 0, ymax = 0, zmax = 0;
@@ -32,7 +32,7 @@ double GeoShape::volume () const
 
   uint32_t y = 2463534242; // seed for random number generation
   int icount = 0; // counter of inside points
-  for (auto i = 0; i < npoints; ++i)
+  for (int i = 0; i < np; ++i)
   {
     // generate three random numbers
     uint32_t x = y;
diff --git a/GeoModelCore/GeoModelKernel/src/GeoShapeIntersection.cxx b/GeoModelCore/GeoModelKernel/src/GeoShapeIntersection.cxx
index 711d2e8c7..219f0109c 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoShapeIntersection.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoShapeIntersection.cxx
@@ -15,10 +15,14 @@ GeoShapeIntersection::GeoShapeIntersection (const GeoShape* A,
                                             const GeoShape* B) :
   m_opA (A) , m_opB (B) {}
 
-
-
-double GeoShapeIntersection::volume () const {
-  return fVolume < 0. ? fVolume = GeoShape::volume() : fVolume.load();
+double GeoShapeIntersection::volume (int npoints) const {
+  double vol = getVolumeValue();
+  if (vol < 0)
+  {
+    vol = GeoShape::volume(npoints);
+    setVolumeValue(vol);
+  }
+  return vol;
 }
 
 void GeoShapeIntersection::extent (double& xmin, double& ymin, double& zmin,
@@ -46,27 +50,26 @@ void GeoShapeIntersection::exec (GeoShapeAction *action) const
   action->getPath ()->push (this);
   action->handleIntersection (this);
   if (action->shouldTerminate ())
-    {
-      action->getPath ()->pop ();
-      return;
-    }
+  {
+    action->getPath ()->pop ();
+    return;
+  }
 
   if (action->getDepthLimit ().isValid ()
       && action->getPath ()->getLength () > action->getDepthLimit ())
-    {
-    }
-
+  {
+  }
   else
   {
     m_opA->exec(action);
     if (action->shouldTerminate ())
-  {
+    {
       action->getPath ()->pop ();
       return;
-  }
+    }
     m_opB->exec(action);
     if (action->shouldTerminate ())
-  {
+    {
       action->getPath ()->pop ();
       return;
     }
diff --git a/GeoModelCore/GeoModelKernel/src/GeoShapeShift.cxx b/GeoModelCore/GeoModelKernel/src/GeoShapeShift.cxx
index 1a6072ca8..373a804a1 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoShapeShift.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoShapeShift.cxx
@@ -64,24 +64,23 @@ void GeoShapeShift::exec (GeoShapeAction *action) const {
   action->getPath ()->push (this);
   action->handleShift (this);
   if (action->shouldTerminate ())
-    {
-      action->getPath ()->pop ();
-      return;
-    }
+  {
+    action->getPath ()->pop ();
+    return;
+  }
 
   if (action->getDepthLimit ().isValid ()
       && action->getPath ()->getLength () > action->getDepthLimit ())
-    {
-    }
-
+  {
+  }
   else
-        {
+  {
     m_op->exec(action);
     if (action->shouldTerminate ())
-        {
+    {
       action->getPath ()->pop ();
       return;
-        }
+    }
   }
   action->getPath()->pop();
 }
diff --git a/GeoModelCore/GeoModelKernel/src/GeoShapeSubtraction.cxx b/GeoModelCore/GeoModelKernel/src/GeoShapeSubtraction.cxx
index 10a919f67..7eebbfbf2 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoShapeSubtraction.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoShapeSubtraction.cxx
@@ -1,8 +1,9 @@
 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "GeoModelKernel/GeoShapeSubtraction.h"
+#include "GeoModelKernel/GeoShapeIntersection.h"
 #include "GeoModelKernel/GeoShapeAction.h"
 #include "GeoModelKernel/GeoPolyhedrizeAction.h"
 #include "GeoModelKernel/GeoPolyhedron.h"
@@ -11,12 +12,35 @@
 const std::string GeoShapeSubtraction::s_classType = "Subtraction";
 const ShapeType GeoShapeSubtraction::s_classTypeID = 0x02;
 
-GeoShapeSubtraction::GeoShapeSubtraction (const GeoShape* A, const GeoShape* B): 
+GeoShapeSubtraction::GeoShapeSubtraction (const GeoShape* A, const GeoShape* B):
     m_opA {A}, m_opB {B} {}
 
-
-double GeoShapeSubtraction::volume () const {
-  return (fVolume < 0.) ? (fVolume = GeoShape::volume()) : fVolume.load();
+double GeoShapeSubtraction::volume (int npoints) const {
+  double vol = getVolumeValue();
+  if (vol < 0)
+  {
+    double xminA, yminA, zminA, xmaxA, ymaxA, zmaxA;
+    double xminB, yminB, zminB, xmaxB, ymaxB, zmaxB;
+    m_opA->extent(xminA, yminA, zminA, xmaxA, ymaxA, zmaxA);
+    m_opB->extent(xminB, yminB, zminB, xmaxB, ymaxB, zmaxB);
+    bool noIntersection =
+      xminA >= xmaxB || yminA >= ymaxB || zminA >= zmaxB ||
+      xminB >= xmaxA || yminB >= ymaxA || zminB >= zmaxA;
+    if (noIntersection) {
+      vol = m_opA->volume(npoints);
+    } else {
+      if (getNoConstituents() > 10) {
+        vol = GeoShape::volume(npoints);
+      } else {
+        GeoIntrusivePtr<GeoShapeIntersection> tmp(new GeoShapeIntersection(m_opA, m_opB));
+        double volumeA = m_opA->volume(npoints);
+        vol = volumeA - tmp->volume(npoints);
+        if (vol < 0.01*volumeA) vol = GeoShape::volume(npoints);
+      }
+    }
+    setVolumeValue(vol);
+  }
+  return vol;
 }
 
 void GeoShapeSubtraction::extent (double& xmin, double& ymin, double& zmin,
@@ -30,34 +54,31 @@ bool GeoShapeSubtraction::contains (double x, double y, double z) const
   return (!getOpA()->contains(x, y, z)) ? false : !getOpB()->contains(x, y, z);
 }
 
-
-
 void GeoShapeSubtraction::exec (GeoShapeAction *action) const
 {
   action->getPath ()->push (this);
   action->handleSubtraction (this);
   if (action->shouldTerminate ())
-    {
-      action->getPath ()->pop ();
-      return;
-    }
+  {
+    action->getPath ()->pop ();
+    return;
+  }
 
   if (action->getDepthLimit ().isValid ()
       && action->getPath ()->getLength () > action->getDepthLimit ())
-    {
-    }
-
+  {
+  }
   else
   {
     m_opA->exec(action);
     if (action->shouldTerminate ())
-  {
+    {
       action->getPath ()->pop ();
       return;
-  }
+    }
     m_opB->exec(action);
     if (action->shouldTerminate ())
-  {
+    {
       action->getPath ()->pop ();
       return;
     }
diff --git a/GeoModelCore/GeoModelKernel/src/GeoShapeUnion.cxx b/GeoModelCore/GeoModelKernel/src/GeoShapeUnion.cxx
index 2acc7c70c..f26055a03 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoShapeUnion.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoShapeUnion.cxx
@@ -1,8 +1,9 @@
 /*
-  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "GeoModelKernel/GeoShapeUnion.h"
+#include "GeoModelKernel/GeoShapeIntersection.h"
 #include "GeoModelKernel/GeoShapeAction.h"
 #include "GeoModelKernel/GeoPolyhedron.h"
 #include "GeoModelKernel/GeoPolyhedrizeAction.h"
@@ -11,13 +12,34 @@
 const std::string GeoShapeUnion::s_classType = "Union";
 const ShapeType GeoShapeUnion::s_classTypeID = 0x01;
 
-GeoShapeUnion::GeoShapeUnion (const GeoShape* A, const GeoShape* B): 
+GeoShapeUnion::GeoShapeUnion (const GeoShape* A, const GeoShape* B):
    m_opA{A}, m_opB{B} {}
 
-
-double GeoShapeUnion::volume () const
+double GeoShapeUnion::volume (int npoints) const
 {
-  return (fVolume < 0.0) ? (fVolume = GeoShape::volume()) : fVolume.load();
+  double vol = getVolumeValue();
+  if (vol < 0)
+  {
+    double xminA, yminA, zminA, xmaxA, ymaxA, zmaxA;
+    double xminB, yminB, zminB, xmaxB, ymaxB, zmaxB;
+    m_opA->extent(xminA, yminA, zminA, xmaxA, ymaxA, zmaxA);
+    m_opB->extent(xminB, yminB, zminB, xmaxB, ymaxB, zmaxB);
+    bool noIntersection =
+      xminA >= xmaxB || yminA >= ymaxB || zminA >= zmaxB ||
+      xminB >= xmaxA || yminB >= ymaxA || zminB >= zmaxA;
+    if (noIntersection) {
+      vol = m_opA->volume(npoints) + m_opB->volume(npoints);
+    } else {
+      if (getNoConstituents() > 10) {
+        vol = GeoShape::volume(npoints);
+      } else {
+        GeoIntrusivePtr<GeoShapeIntersection> tmp(new GeoShapeIntersection(m_opA, m_opB));
+        vol = m_opA->volume(npoints) + m_opB->volume(npoints) - tmp->volume(npoints);
+      }
+    }
+    setVolumeValue(vol);
+  }
+  return vol;
 }
 
 void GeoShapeUnion::extent (double& xmin, double& ymin, double& zmin,
@@ -44,27 +66,26 @@ void GeoShapeUnion::exec (GeoShapeAction *action) const
   action->getPath ()->push (this);
   action->handleUnion (this);
   if (action->shouldTerminate ())
-    {
-      action->getPath ()->pop ();
-      return;
-    }
+  {
+    action->getPath ()->pop ();
+    return;
+  }
 
   if (action->getDepthLimit ().isValid ()
       && action->getPath ()->getLength () > action->getDepthLimit ())
-    {
-    }
-
+  {
+  }
   else
   {
     m_opA->exec(action);
     if (action->shouldTerminate ())
-  {
+    {
       action->getPath ()->pop ();
       return;
-  }
+    }
     m_opB->exec(action);
     if (action->shouldTerminate ())
-  {
+    {
       action->getPath ()->pop ();
       return;
     }
diff --git a/GeoModelCore/GeoModelKernel/src/GeoSimplePolygonBrep.cxx b/GeoModelCore/GeoModelKernel/src/GeoSimplePolygonBrep.cxx
index 2e9eb8de4..001088ace 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoSimplePolygonBrep.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoSimplePolygonBrep.cxx
@@ -12,11 +12,10 @@ const std::string GeoSimplePolygonBrep::s_classType = "SimplePolygonBrep";
 const ShapeType GeoSimplePolygonBrep::s_classTypeID = 0x20;
 
 GeoSimplePolygonBrep::GeoSimplePolygonBrep(double dz)
-   : m_dZ(dz){}
+   : m_dZ(dz)
+{}
 
-
-
-double GeoSimplePolygonBrep::volume () const {
+double GeoSimplePolygonBrep::volume (int) const {
   if (!isValid())
     throw std::runtime_error ("Volume requested for incomplete simple polygon brep");
   int n = getNVertices();
@@ -64,7 +63,6 @@ bool GeoSimplePolygonBrep::contains (double x, double y, double z) const
   return in;
 }
 
-
 void GeoSimplePolygonBrep::addVertex(double XVertex, double YVertex)
 {
   m_xVertices.push_back(XVertex);
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTessellatedSolid.cxx b/GeoModelCore/GeoModelKernel/src/GeoTessellatedSolid.cxx
index 2f7a3fabc..19fcfee86 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoTessellatedSolid.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoTessellatedSolid.cxx
@@ -13,7 +13,7 @@ GeoTessellatedSolid::GeoTessellatedSolid()
 {
 }
 
-double GeoTessellatedSolid::volume() const
+double GeoTessellatedSolid::volume(int) const
 {
   if (!isValid ())
     throw std::runtime_error ("Volume requested for incomplete tessellated solid");
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTorus.cxx b/GeoModelCore/GeoModelKernel/src/GeoTorus.cxx
index 717d58083..8626a6cb1 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoTorus.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoTorus.cxx
@@ -15,11 +15,9 @@ GeoTorus::GeoTorus (double Rmin, double Rmax, double Rtor, double SPhi, double D
    , m_rTor (Rtor)
    , m_sPhi (SPhi)
    , m_dPhi (DPhi)
-{
-}
-
+{}
 
-double GeoTorus::volume () const
+double GeoTorus::volume (int) const
 {
 #ifndef M_PI
   constexpr double M_PI = 3.14159265358979323846;
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTrap.cxx b/GeoModelCore/GeoModelKernel/src/GeoTrap.cxx
index 99918eff3..45b39da98 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoTrap.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoTrap.cxx
@@ -20,11 +20,9 @@ GeoTrap::GeoTrap (double ZHalfLength, double Theta, double Phi, double Dydzn, do
   , m_dxdyndzp (Dxdyndzp)
   , m_dxdypdzp (Dxdypdzp)
   , m_angleydzp (Angleydzp)
-{
-}
+{}
 
-
-double GeoTrap::volume () const
+double GeoTrap::volume (int) const
 {
   double dz = m_zHalfLength;
   double dy1 = m_dydzn;
@@ -108,7 +106,6 @@ bool GeoTrap::contains (double x, double y, double z) const
   return (std::abs(x0) - dx <= 0.0);
 }
 
-
 void GeoTrap::exec (GeoShapeAction *action) const
 {
   action->handleTrap(this);
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTrd.cxx b/GeoModelCore/GeoModelKernel/src/GeoTrd.cxx
index 1e0b40c03..3d35073e9 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoTrd.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoTrd.cxx
@@ -16,8 +16,7 @@ GeoTrd::GeoTrd (double XHalfLength1, double XHalfLength2, double YHalfLength1, d
   , m_zHalfLength (ZHalfLength)
 {}
 
-
-double GeoTrd::volume () const
+double GeoTrd::volume (int) const
 {
   double fDz = m_zHalfLength;
   double fDy1 = m_yHalfLength1;
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTube.cxx b/GeoModelCore/GeoModelKernel/src/GeoTube.cxx
index b05760ed7..35f10c2fa 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoTube.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoTube.cxx
@@ -14,14 +14,9 @@ GeoTube::GeoTube (double RMin, double RMax, double ZHalfLength)
   : m_rMin (RMin)
   , m_rMax (RMax)
   , m_zHalfLength (ZHalfLength)
-{
-}
-
-GeoTube::~GeoTube()
-{
-}
+{}
 
-double GeoTube::volume () const
+double GeoTube::volume (int) const
 {
 #ifndef M_PI
   constexpr double M_PI = 3.14159265358979323846;
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTubs.cxx b/GeoModelCore/GeoModelKernel/src/GeoTubs.cxx
index 0a7776b96..ea699ee07 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoTubs.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoTubs.cxx
@@ -15,11 +15,9 @@ GeoTubs::GeoTubs (double RMin, double RMax, double ZHalfLength, double SPhi, dou
   , m_zHalfLength (ZHalfLength)
   , m_sPhi (SPhi)
   , m_dPhi (DPhi)
-{
-}
+{}
 
-
-double GeoTubs::volume () const
+double GeoTubs::volume (int) const
 {
   return m_dPhi * (m_rMax * m_rMax - m_rMin * m_rMin) * m_zHalfLength;
 }
@@ -50,7 +48,6 @@ bool GeoTubs::contains (double x, double y, double z) const
   return (m_dPhi <= M_PI) ? (ds <= 0 && de <= 0) : (ds <= 0 || de <= 0);
 }
 
-
 void GeoTubs::exec (GeoShapeAction *action) const
 {
   action->handleTubs(this);
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTwistedTrap.cxx b/GeoModelCore/GeoModelKernel/src/GeoTwistedTrap.cxx
index 07e3c994b..35311e4d9 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoTwistedTrap.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoTwistedTrap.cxx
@@ -85,7 +85,7 @@ GeoTwistedTrap(double  pPhiTwist,  // twist angle
 }
 
 
-double GeoTwistedTrap::volume () const
+double GeoTwistedTrap::volume (int) const
 {
   return m_dz * ((m_dx1 + m_dx2 + m_dx3 + m_dx4) * (m_dy1 + m_dy2) +
                  (m_dx4 + m_dx3 - m_dx2 - m_dx1) * (m_dy2 - m_dy1) * (1./3.));
-- 
GitLab