Repository: NREL/OpenStudio (language: python)

Issue:
Fix volume calculation to match EnergyPlus (Bugzilla #730)
On 2012-06-08 14:19:09, @DavidGoldwasser wrote:

> Just wanted to file this before I forgot about it. Current method will have problem with spaces that have sloped walls or ceilings.

On 2012-09-09 23:24:43, @DavidGoldwasser wrote:

> I added BadSimResults but that would be an edge case where walls or roofs are not vertical.

On 2012-12-17 13:42:44, @DavidGoldwasser wrote:

> This was inadvertently under C#

Keywords: BadSimResults

Fix volume calculation to match EnergyPlus (Bugzilla #730)
On 2012-06-08 14:19:09, @DavidGoldwasser wrote:

> Just wanted to file this before I forgot about it. Current method will have problem with spaces that have sloped walls or ceilings.

On 2012-09-09 23:24:43, @DavidGoldwasser wrote:

> I added BadSimResults but that would be an edge case where walls or roofs are not vertical.

On 2012-12-17 13:42:44, @DavidGoldwasser wrote:

> This was inadvertently under C#

Keywords: BadSimResults

Fix volume calculation to match EnergyPlus (Bugzilla #730)
On 2012-06-08 14:19:09, @DavidGoldwasser wrote:

> Just wanted to file this before I forgot about it. Current method will have problem with spaces that have sloped walls or ceilings.

On 2012-09-09 23:24:43, @DavidGoldwasser wrote:

> I added BadSimResults but that would be an edge case where walls or roofs are not vertical.

On 2012-12-17 13:42:44, @DavidGoldwasser wrote:

> This was inadvertently under C#

Keywords: BadSimResults



Discussion / hints:
@macumber I'll keep this open, but I know it is more of a feature than a bug, but it can be important for proper calculation of some internal load values on forward translation of complex geometry.

I propose using the divergence theorem to compute the space volume, instead of multiplying the space area by a "mean" space height.
![image](https://user-images.githubusercontent.com/38532229/39042198-89e2f86e-448a-11e8-9afb-04f47bf9f674.png)
Then, using that every surface belongs to a plane, whose normal is parallel to its coefficients a, b and c;
![image](https://user-images.githubusercontent.com/38532229/39042288-b9944770-448a-11e8-9fa3-f9782df3947a.png)
only a constant integral is needed in every surface of the space.
![image](https://user-images.githubusercontent.com/38532229/39042362-e27e2034-448a-11e8-8453-52c764eab314.png)
Then, the volume integral reduces to
![image](https://user-images.githubusercontent.com/38532229/39042456-2891eb00-448b-11e8-8aaf-37661a4adaf0.png)
where d_i is the d coefficent of the plane where the i-th surface belongs and A_i is the area of the i-th surface.

I've tested it using this ruby code.
![image](https://user-images.githubusercontent.com/38532229/39042540-5e6ecc34-448b-11e8-84aa-fb9d3adf4804.png)


 

@agonzalezesteve thank you, this is awesome. 

I just tested the following single space geometry: 
![image](https://user-images.githubusercontent.com/10799524/39071790-9a1af36a-44a5-11e8-9c28-0fb74ff1fad2.png)

By hand, I calculate the volume at 17.6642 m^3 (measuring from sketchup in ft/in and converting ft^3 to m^3 in excel)
The default `space.volume` method for this space returns 19.0972 m^3
Your method returns 17.5825 m^3.
I'd say that's an improvement!

I've had many problems with OpenStudio geometry from sketchup ending up with EnergyPlus 'volume not enclosed' errors, to the point where I've taken to hard-assigning zone volume from a script. For complex spaces, I've had to hand-calc the volumes. I think this method is a huge improvement, and the dev team might consider implementing it to set zone volume at forward translation to avoid problems that EnergyPlus has calculating non-enclosed zone volumes with tiny surface gaps. 


The divergence calculation does not work for volumes which are not fully enclosed, I didn't implement it for OpenStudio because we don't have a check that the volumes are fully enclosed (no gaps or overlapping surfaces).  I agree that it is best to match E+, I have been holding out for a solution where OS and E+ *share* code for calculations like this.  That way we can make sure we are using the best available algorithms.
@macumber ok maybe not implement it in OS by default, but the above algorithm is still better than zone volumes being set to 10 m^3 when E+ detects an edge that's not shared by two surfaces due to tiny gaps from sketchup. It will "work", just not "totally correctly" 😄 
I'm not really following @eringold, it sounds like you are having problems with the EnergyPlus calculation which fails due to tiny gaps (which invalidate the divergence calculation).  Is that right? Is your problem with the OS calculation or the E+ one?
@macumber both.  I've had problems with EnergyPlus failing to calculate volume from OS geometry that I just could not get  to be fully enclosed when translated*, no matter how careful I was drawing it in OpenStudio. Hence the desire to hard-size zone volumes. But the easy way to do that is with the OpenStudio volumes, which can fail for funny-shaped zones like the one above (though I think the EnergyPlus volume calc might be more robust for that? - if it's fully enclosed). So I like the divergence method because, even though it's not perfect for unenclosed volumes (with, like, tiny gaps), the script above will still run and return a value close to the 'actual' volume, but can also handle complex geometry. 

(I know it would be **_best_** to not have complex geometry in models to begin with, but stuff happens..)

*I meant to look into this in more detail but had a timeline that wouldn't allow for it.
@macumber in which cases the geometry is not fully enclosed?

It is true that the divergence calculation doesn't work with unenclosed spaces, but the actual method neither if there is no floor or roof.

Furthermore, the actual method is a product of a "mean" height times the floor area, that in complex geometries this simplification fails.

In short, the current and the proposed methods for volume computation don't work with unenclosed volumes but the divergence method is more accurate in complex geometries.
First of all, thank you @agonzalezesteve for this conversation, it is great to have you contributing here!

I don't know all the cases that geometry isn't fully enclosed.  Thank you very much for your suggestion, I agree with your evaluation of the two approaches, the only thing I am unsure about is the robustness of both approaches.  I was under the impression that you could get wildy wrong volumes from the divergence routine if you were missing a wall or something like that so the current routine is more robust to failure, I could be wrong here.  In the end I suspect that a hybrid routine (e.g. use divergence if fully enclosed, issue warning and use some other method if not) would be best suited here.  My hope was to implement this hybrid routine in a shared library used by both OpenStudio and EnergyPlus, I still think this is possible.

Currently OpenStudio does not hard set volume on translation to E+ so I don't think OpenStudio's current volume method affects most users.  Am I missing some case where OpenStudio's volume calculation is commonly invoked?  If a user wants to hard set volume in a measure then I think it is good for them to be able to choose between methods.  I don't have a strong opinion on which one should be the default so my preference is to keep the current algorithm for now and point users to your snippet if they want to try that approach.
@macumber thank you for the welcome!

You are right, the divergence method fails if one wall is missing and the current method is more robust in these cases.

I just tried to find a new volume algorithm because the current method gave me wrong values of air change per hour in spaces with complex geometry, like staircases inside buildings.
https://github.com/NREL/OpenStudio/blob/a007b7801e3bb05731f012f585facdd814ef31fb/src/model/Space.cpp#L887-L916
The corresponding entry point for E+ calculation is around here: https://github.com/NREL/EnergyPlus/blob/29b5decc089394689a97932f66fa64d5a2fac062/src/EnergyPlus/SurfaceGeometry.cc#L11506-L11534
This has not been fixed since 2013! I tried the workaround osThermalZone.setVolume but, while this fixes it for the .OSM file, it does not fix it for the .gbxml output which continues to calculate the wrong volume for more complex shapes (mainly slanted hip roofs etc).
@macumber I'll keep this open, but I know it is more of a feature than a bug, but it can be important for proper calculation of some internal load values on forward translation of complex geometry.

I propose using the divergence theorem to compute the space volume, instead of multiplying the space area by a "mean" space height.
![image](https://user-images.githubusercontent.com/38532229/39042198-89e2f86e-448a-11e8-9afb-04f47bf9f674.png)
Then, using that every surface belongs to a plane, whose normal is parallel to its coefficients a, b and c;
![image](https://user-images.githubusercontent.com/38532229/39042288-b9944770-448a-11e8-9fa3-f9782df3947a.png)
only a constant integral is needed in every surface of the space.
![image](https://user-images.githubusercontent.com/38532229/39042362-e27e2034-448a-11e8-8453-52c764eab314.png)
Then, the volume integral reduces to
![image](https://user-images.githubusercontent.com/38532229/39042456-2891eb00-448b-11e8-8aaf-37661a4adaf0.png)
where d_i is the d coefficent of the plane where the i-th surface belongs and A_i is the area of the i-th surface.

I've tested it using this ruby code.
![image](https://user-images.githubusercontent.com/38532229/39042540-5e6ecc34-448b-11e8-84aa-fb9d3adf4804.png)


 

@agonzalezesteve thank you, this is awesome. 

I just tested the following single space geometry: 
![image](https://user-images.githubusercontent.com/10799524/39071790-9a1af36a-44a5-11e8-9c28-0fb74ff1fad2.png)

By hand, I calculate the volume at 17.6642 m^3 (measuring from sketchup in ft/in and converting ft^3 to m^3 in excel)
The default `space.volume` method for this space returns 19.0972 m^3
Your method returns 17.5825 m^3.
I'd say that's an improvement!

I've had many problems with OpenStudio geometry from sketchup ending up with EnergyPlus 'volume not enclosed' errors, to the point where I've taken to hard-assigning zone volume from a script. For complex spaces, I've had to hand-calc the volumes. I think this method is a huge improvement, and the dev team might consider implementing it to set zone volume at forward translation to avoid problems that EnergyPlus has calculating non-enclosed zone volumes with tiny surface gaps. 


The divergence calculation does not work for volumes which are not fully enclosed, I didn't implement it for OpenStudio because we don't have a check that the volumes are fully enclosed (no gaps or overlapping surfaces).  I agree that it is best to match E+, I have been holding out for a solution where OS and E+ *share* code for calculations like this.  That way we can make sure we are using the best available algorithms.
@macumber ok maybe not implement it in OS by default, but the above algorithm is still better than zone volumes being set to 10 m^3 when E+ detects an edge that's not shared by two surfaces due to tiny gaps from sketchup. It will "work", just not "totally correctly" 😄 
I'm not really following @eringold, it sounds like you are having problems with the EnergyPlus calculation which fails due to tiny gaps (which invalidate the divergence calculation).  Is that right? Is your problem with the OS calculation or the E+ one?
@macumber both.  I've had problems with EnergyPlus failing to calculate volume from OS geometry that I just could not get  to be fully enclosed when translated*, no matter how careful I was drawing it in OpenStudio. Hence the desire to hard-size zone volumes. But the easy way to do that is with the OpenStudio volumes, which can fail for funny-shaped zones like the one above (though I think the EnergyPlus volume calc might be more robust for that? - if it's fully enclosed). So I like the divergence method because, even though it's not perfect for unenclosed volumes (with, like, tiny gaps), the script above will still run and return a value close to the 'actual' volume, but can also handle complex geometry. 

(I know it would be **_best_** to not have complex geometry in models to begin with, but stuff happens..)

*I meant to look into this in more detail but had a timeline that wouldn't allow for it.
@macumber in which cases the geometry is not fully enclosed?

It is true that the divergence calculation doesn't work with unenclosed spaces, but the actual method neither if there is no floor or roof.

Furthermore, the actual method is a product of a "mean" height times the floor area, that in complex geometries this simplification fails.

In short, the current and the proposed methods for volume computation don't work with unenclosed volumes but the divergence method is more accurate in complex geometries.
First of all, thank you @agonzalezesteve for this conversation, it is great to have you contributing here!

I don't know all the cases that geometry isn't fully enclosed.  Thank you very much for your suggestion, I agree with your evaluation of the two approaches, the only thing I am unsure about is the robustness of both approaches.  I was under the impression that you could get wildy wrong volumes from the divergence routine if you were missing a wall or something like that so the current routine is more robust to failure, I could be wrong here.  In the end I suspect that a hybrid routine (e.g. use divergence if fully enclosed, issue warning and use some other method if not) would be best suited here.  My hope was to implement this hybrid routine in a shared library used by both OpenStudio and EnergyPlus, I still think this is possible.

Currently OpenStudio does not hard set volume on translation to E+ so I don't think OpenStudio's current volume method affects most users.  Am I missing some case where OpenStudio's volume calculation is commonly invoked?  If a user wants to hard set volume in a measure then I think it is good for them to be able to choose between methods.  I don't have a strong opinion on which one should be the default so my preference is to keep the current algorithm for now and point users to your snippet if they want to try that approach.
@macumber thank you for the welcome!

You are right, the divergence method fails if one wall is missing and the current method is more robust in these cases.

I just tried to find a new volume algorithm because the current method gave me wrong values of air change per hour in spaces with complex geometry, like staircases inside buildings.
https://github.com/NREL/OpenStudio/blob/a007b7801e3bb05731f012f585facdd814ef31fb/src/model/Space.cpp#L887-L916
The corresponding entry point for E+ calculation is around here: https://github.com/NREL/EnergyPlus/blob/29b5decc089394689a97932f66fa64d5a2fac062/src/EnergyPlus/SurfaceGeometry.cc#L11506-L11534
This has not been fixed since 2013! I tried the workaround osThermalZone.setVolume but, while this fixes it for the .OSM file, it does not fix it for the .gbxml output which continues to calculate the wrong volume for more complex shapes (mainly slanted hip roofs etc).
@macumber I'll keep this open, but I know it is more of a feature than a bug, but it can be important for proper calculation of some internal load values on forward translation of complex geometry.

I propose using the divergence theorem to compute the space volume, instead of multiplying the space area by a "mean" space height.
![image](https://user-images.githubusercontent.com/38532229/39042198-89e2f86e-448a-11e8-9afb-04f47bf9f674.png)
Then, using that every surface belongs to a plane, whose normal is parallel to its coefficients a, b and c;
![image](https://user-images.githubusercontent.com/38532229/39042288-b9944770-448a-11e8-9fa3-f9782df3947a.png)
only a constant integral is needed in every surface of the space.
![image](https://user-images.githubusercontent.com/38532229/39042362-e27e2034-448a-11e8-8453-52c764eab314.png)
Then, the volume integral reduces to
![image](https://user-images.githubusercontent.com/38532229/39042456-2891eb00-448b-11e8-8aaf-37661a4adaf0.png)
where d_i is the d coefficent of the plane where the i-th surface belongs and A_i is the area of the i-th surface.

I've tested it using this ruby code.
![image](https://user-images.githubusercontent.com/38532229/39042540-5e6ecc34-448b-11e8-84aa-fb9d3adf4804.png)


 

@agonzalezesteve thank you, this is awesome. 

I just tested the following single space geometry: 
![image](https://user-images.githubusercontent.com/10799524/39071790-9a1af36a-44a5-11e8-9c28-0fb74ff1fad2.png)

By hand, I calculate the volume at 17.6642 m^3 (measuring from sketchup in ft/in and converting ft^3 to m^3 in excel)
The default `space.volume` method for this space returns 19.0972 m^3
Your method returns 17.5825 m^3.
I'd say that's an improvement!

I've had many problems with OpenStudio geometry from sketchup ending up with EnergyPlus 'volume not enclosed' errors, to the point where I've taken to hard-assigning zone volume from a script. For complex spaces, I've had to hand-calc the volumes. I think this method is a huge improvement, and the dev team might consider implementing it to set zone volume at forward translation to avoid problems that EnergyPlus has calculating non-enclosed zone volumes with tiny surface gaps. 


The divergence calculation does not work for volumes which are not fully enclosed, I didn't implement it for OpenStudio because we don't have a check that the volumes are fully enclosed (no gaps or overlapping surfaces).  I agree that it is best to match E+, I have been holding out for a solution where OS and E+ *share* code for calculations like this.  That way we can make sure we are using the best available algorithms.
@macumber ok maybe not implement it in OS by default, but the above algorithm is still better than zone volumes being set to 10 m^3 when E+ detects an edge that's not shared by two surfaces due to tiny gaps from sketchup. It will "work", just not "totally correctly" 😄 
I'm not really following @eringold, it sounds like you are having problems with the EnergyPlus calculation which fails due to tiny gaps (which invalidate the divergence calculation).  Is that right? Is your problem with the OS calculation or the E+ one?
@macumber both.  I've had problems with EnergyPlus failing to calculate volume from OS geometry that I just could not get  to be fully enclosed when translated*, no matter how careful I was drawing it in OpenStudio. Hence the desire to hard-size zone volumes. But the easy way to do that is with the OpenStudio volumes, which can fail for funny-shaped zones like the one above (though I think the EnergyPlus volume calc might be more robust for that? - if it's fully enclosed). So I like the divergence method because, even though it's not perfect for unenclosed volumes (with, like, tiny gaps), the script above will still run and return a value close to the 'actual' volume, but can also handle complex geometry. 

(I know it would be **_best_** to not have complex geometry in models to begin with, but stuff happens..)

*I meant to look into this in more detail but had a timeline that wouldn't allow for it.
@macumber in which cases the geometry is not fully enclosed?

It is true that the divergence calculation doesn't work with unenclosed spaces, but the actual method neither if there is no floor or roof.

Furthermore, the actual method is a product of a "mean" height times the floor area, that in complex geometries this simplification fails.

In short, the current and the proposed methods for volume computation don't work with unenclosed volumes but the divergence method is more accurate in complex geometries.
First of all, thank you @agonzalezesteve for this conversation, it is great to have you contributing here!

I don't know all the cases that geometry isn't fully enclosed.  Thank you very much for your suggestion, I agree with your evaluation of the two approaches, the only thing I am unsure about is the robustness of both approaches.  I was under the impression that you could get wildy wrong volumes from the divergence routine if you were missing a wall or something like that so the current routine is more robust to failure, I could be wrong here.  In the end I suspect that a hybrid routine (e.g. use divergence if fully enclosed, issue warning and use some other method if not) would be best suited here.  My hope was to implement this hybrid routine in a shared library used by both OpenStudio and EnergyPlus, I still think this is possible.

Currently OpenStudio does not hard set volume on translation to E+ so I don't think OpenStudio's current volume method affects most users.  Am I missing some case where OpenStudio's volume calculation is commonly invoked?  If a user wants to hard set volume in a measure then I think it is good for them to be able to choose between methods.  I don't have a strong opinion on which one should be the default so my preference is to keep the current algorithm for now and point users to your snippet if they want to try that approach.
@macumber thank you for the welcome!

You are right, the divergence method fails if one wall is missing and the current method is more robust in these cases.

I just tried to find a new volume algorithm because the current method gave me wrong values of air change per hour in spaces with complex geometry, like staircases inside buildings.
https://github.com/NREL/OpenStudio/blob/a007b7801e3bb05731f012f585facdd814ef31fb/src/model/Space.cpp#L887-L916
The corresponding entry point for E+ calculation is around here: https://github.com/NREL/EnergyPlus/blob/29b5decc089394689a97932f66fa64d5a2fac062/src/EnergyPlus/SurfaceGeometry.cc#L11506-L11534
This has not been fixed since 2013! I tried the workaround osThermalZone.setVolume but, while this fixes it for the .OSM file, it does not fix it for the .gbxml output which continues to calculate the wrong volume for more complex shapes (mainly slanted hip roofs etc).

Fix patch:
diff --git a/src/model/ModelGeometry.i b/src/model/ModelGeometry.i
index b094d68c133..3ae7600afb7 100644
--- a/src/model/ModelGeometry.i
+++ b/src/model/ModelGeometry.i
@@ -14,6 +14,7 @@
 %{
   #include <utilities/geometry/Transformation.hpp>
   #include <utilities/geometry/BoundingBox.hpp>
+  #include <utilities/geometry/Polyhedron.hpp>
   #include <utilities/data/TimeSeries.hpp>
   #include <utilities/sql/SqlFile.hpp>
 
diff --git a/src/model/Space.cpp b/src/model/Space.cpp
index 853802f2865..4789338e1fa 100644
--- a/src/model/Space.cpp
+++ b/src/model/Space.cpp
@@ -127,6 +127,7 @@
 #include "../utilities/geometry/EulerAngles.hpp"
 #include "../utilities/geometry/BoundingBox.hpp"
 #include "../utilities/geometry/Polygon3d.hpp"
+#include "../utilities/geometry/Polyhedron.hpp"
 
 #include "../utilities/core/Assert.hpp"
 
@@ -901,12 +902,42 @@ namespace model {
       return result;
     }
 
+    Polyhedron Space_Impl::polyhedron() const {
+      std::vector<Surface3d> surface3ds;
+      for (auto& surface : surfaces()) {
+        surface3ds.emplace_back(surface.vertices(), surface.nameString());
+      }
+      return {surface3ds};
+    }
+
+    bool Space_Impl::isEnclosedVolume() const {
+      auto volumePoly = this->polyhedron();
+      auto [isVolEnclosed, edgesNot2] = volumePoly.isEnclosedVolume();
+      if (!isVolEnclosed) {
+        LOG(Warn, briefDescription() << " is not enclosed, there are " << edgesNot2.size() << " edges that aren't used exactly twice");
+        for (const Surface3dEdge& edge : edgesNot2) {
+          LOG(Debug, edge);
+        }
+      }
+      return isVolEnclosed;
+    }
+
     double Space_Impl::volume() const {
       boost::optional<double> value = getDouble(OS_SpaceFields::Volume, true);
       if (value) {
         return value.get();
       }
 
+      auto volumePoly = this->polyhedron();
+
+      auto [isVolEnclosed, edgesNot2] = volumePoly.isEnclosedVolume();
+      if (isVolEnclosed) {
+        return volumePoly.calcPolyhedronVolume();
+      }
+
+      LOG(Warn, briefDescription() << " is not enclosed, there are " << edgesNot2.size()
+                                   << " edges that aren't used exactly twice. Volume calculation will be potentially inaccurate");
+
       double result = 0;
 
       // TODO: need a better method
@@ -3334,6 +3365,14 @@ namespace model {
     return getImpl<detail::Space_Impl>()->exposedPerimeter(buildingPerimeter);
   }
 
+  Polyhedron Space::polyhedron() const {
+    return getImpl<detail::Space_Impl>()->polyhedron();
+  }
+
+  bool Space::isEnclosedVolume() const {
+    return getImpl<detail::Space_Impl>()->isEnclosedVolume();
+  }
+
   /// @cond
   Space::Space(std::shared_ptr<detail::Space_Impl> impl) : PlanarSurfaceGroup(std::move(impl)) {}
   /// @endcond
diff --git a/src/model/Space.hpp b/src/model/Space.hpp
index 9f024f5c7f5..351668e85e9 100644
--- a/src/model/Space.hpp
+++ b/src/model/Space.hpp
@@ -34,7 +34,9 @@
 #include "PlanarSurfaceGroup.hpp"
 
 namespace openstudio {
+
 class Polygon3d;
+class Polyhedron;
 
 namespace model {
 
@@ -596,6 +598,9 @@ namespace model {
     // Calculates the exposed perimeter of a space
     double exposedPerimeter(const Polygon3d& buildingPerimeter) const;
 
+    Polyhedron polyhedron() const;
+    bool isEnclosedVolume() const;
+
     //@}
    protected:
     /// @cond
diff --git a/src/model/Space_Impl.hpp b/src/model/Space_Impl.hpp
index cd9f0ab4d90..289f9977470 100644
--- a/src/model/Space_Impl.hpp
+++ b/src/model/Space_Impl.hpp
@@ -37,6 +37,10 @@
 #include <boost/geometry/geometries/adapted/boost_tuple.hpp>
 
 namespace openstudio {
+
+class Polygon3d;
+class Polyhedron;
+
 namespace model {
 
   // forward declarations
@@ -478,6 +482,9 @@ namespace model {
 
       double exposedPerimeter(const Polygon3d& buildingPerimeter) const;
 
+      Polyhedron polyhedron() const;
+      bool isEnclosedVolume() const;
+
      private:
       REGISTER_LOGGER("openstudio.model.Space");
 
diff --git a/src/utilities/CMakeLists.txt b/src/utilities/CMakeLists.txt
index 4df09af97ae..99935d60f78 100644
--- a/src/utilities/CMakeLists.txt
+++ b/src/utilities/CMakeLists.txt
@@ -143,6 +143,8 @@ set(geometry_src
   geometry/Vector3d.cpp
   geometry/Polygon3d.hpp
   geometry/Polygon3d.cpp
+  geometry/Polyhedron.hpp
+  geometry/Polyhedron.cpp
   ../polypartition/polypartition.cpp
 )
 
@@ -315,6 +317,7 @@ set(${target_name}_test_src
   geometry/Test/ThreeJS_GTest.cpp
   geometry/Test/FloorplanJS_GTest.cpp
   geometry/Test/Transformation_GTest.cpp
+  geometry/Test/Polyhedron_GTest.cpp
 
   math/test/FloatCompare_GTest.cpp
   math/test/Permutation_GTest.cpp
diff --git a/src/utilities/geometry/Geometry.cpp b/src/utilities/geometry/Geometry.cpp
index d15029a1cff..409e028d4bf 100644
--- a/src/utilities/geometry/Geometry.cpp
+++ b/src/utilities/geometry/Geometry.cpp
@@ -850,4 +850,21 @@ bool applyViewAndDaylightingGlassRatios(double viewGlassToWallRatio, double dayl
   return true;
 }
 
+bool isAlmostEqual3dPt(const Point3d& lhs, const Point3d& rhs, double tol) {
+  // TODO: this is what E+ does... I think I would prefer just calling getDistance...
+  return ((std::abs(lhs.x() - rhs.x()) < tol) && (std::abs(lhs.y() - rhs.y()) < tol) && (std::abs(lhs.z() - rhs.z()) < tol));
+  // return getDistance(lhs, rhs) < tol;
+}
+
+bool isPointOnLineBetweenPoints(const Point3d& start, const Point3d& end, const Point3d& test, double tol) {
+  // The tolerance has to be low enough. Take for eg a plenum that has an edge that's 30meters long, you risk adding point from the floor to
+  // the roof, cf E+ #7383
+  // compute the shortest distance from the point to the line first to avoid false positive
+  double distance = getDistancePointToLineSegment(test, {start, end});
+  if (distance < tol) {  // getDistancePointToLineSegment always positive, it's calculated as norml_L2
+    return (std::abs((getDistance(start, end) - (getDistance(start, test) + getDistance(test, end)))) < tol);
+  }
+  return false;
+}
+
 }  // namespace openstudio
diff --git a/src/utilities/geometry/Geometry.hpp b/src/utilities/geometry/Geometry.hpp
index e816b24fdd5..3a3540279d1 100644
--- a/src/utilities/geometry/Geometry.hpp
+++ b/src/utilities/geometry/Geometry.hpp
@@ -127,6 +127,11 @@ UTILITIES_API bool applyViewAndDaylightingGlassRatios(double viewGlassToWallRati
                                                       std::vector<Point3d>& daylightingVertices, std::vector<Point3d>& exteriorShadingVertices,
                                                       std::vector<Point3d>& interiorShelfVertices);
 
+// Checks that a point is **almost** equal to another (with some tolerance)
+UTILITIES_API bool isAlmostEqual3dPt(const Point3d& lhs, const Point3d& rhs, double tol = 0.0127);
+
+UTILITIES_API bool isPointOnLineBetweenPoints(const Point3d& start, const Point3d& end, const Point3d& test, double tol = 0.0127);
+
 }  // namespace openstudio
 
 #endif  //UTILITIES_GEOMETRY_GEOMETRY_HPP
diff --git a/src/utilities/geometry/Geometry.i b/src/utilities/geometry/Geometry.i
index 951da8269bd..25217dd3d90 100644
--- a/src/utilities/geometry/Geometry.i
+++ b/src/utilities/geometry/Geometry.i
@@ -24,6 +24,7 @@
   #include <utilities/geometry/ThreeJS.hpp>
   #include <utilities/geometry/FloorplanJS.hpp>
   #include <utilities/geometry/RoofGeometry.hpp>
+  #include <utilities/geometry/Polyhedron.hpp>
 
   #include <utilities/units/Quantity.hpp>
   #include <utilities/units/Unit.hpp>
@@ -54,6 +55,7 @@
 %template(OptionalFloorplanJS) boost::optional<openstudio::FloorplanJS>;
 %template(OptionalFloorplanObject) boost::optional<openstudio::FloorplanObject>;
 %template(OptionalPolygon3d) boost::optional<openstudio::Polygon3d>;
+%template(OptionalPolyhedron) boost::optional<openstudio::Polyhedron>;
 
 // create an instantiation of the vector classes
 // Note JM 2019-04-16: No need to ignore std::vector<T>::vector/resize when you have a default constructor
@@ -64,6 +66,8 @@
 %template(PointLatLonVector) std::vector<openstudio::PointLatLon>;
 %template(Vector3dVector) std::vector<openstudio::Vector3d>;
 %template(Polygon3dVector) std::vector<openstudio::Polygon3d>;
+%template(PolyhedronVector) std::vector<openstudio::Polyhedron>;
+%template(Surface3dVector) std::vector<openstudio::Surface3d>;
 
 %ignore std::vector<openstudio::Plane>::vector(size_type);
 %ignore std::vector<openstudio::Plane>::resize(size_type);
@@ -92,6 +96,10 @@
 %ignore std::vector<openstudio::FloorplanObject>::resize(size_type);
 %template(FloorplanObjectVector) std::vector<openstudio::FloorplanObject>;
 
+%ignore std::vector<openstudio::Surface3dEdge>::vector(size_type);
+%ignore std::vector<openstudio::Surface3dEdge>::resize(size_type);
+%template(Surface3dEdgeVector) std::vector<openstudio::Surface3dEdge>;
+
 %template(SizeTVector) std::vector<size_t>;
 %template(StringStringMap) std::map<std::string, std::string>;
 
@@ -110,6 +118,7 @@
 %include <utilities/geometry/ThreeJS.hpp>
 %include <utilities/geometry/FloorplanJS.hpp>
 %include <utilities/geometry/RoofGeometry.hpp>
+%include <utilities/geometry/Polyhedron.hpp>
 
 %extend openstudio::Vector3d{
   std::string __str__() const {
@@ -135,6 +144,22 @@
   }
 }
 
+%extend openstudio::Plane {
+  std::string __str__() const {
+    std::ostringstream os;
+    os << *self;
+    return os.str();
+  }
+}
+
+%extend openstudio::Surface3dEdge {
+  std::string __str__() const {
+    std::ostringstream os;
+    os << *self;
+    return os.str();
+  }
+}
+
 %extend openstudio::Transformation {
 
   std::string __str__() const {
diff --git a/src/utilities/geometry/Polyhedron.cpp b/src/utilities/geometry/Polyhedron.cpp
new file mode 100644
index 00000000000..b4c32121323
--- /dev/null
+++ b/src/utilities/geometry/Polyhedron.cpp
@@ -0,0 +1,291 @@
+/***********************************************************************************************************************
+*  OpenStudio(R), Copyright (c) 2008-2022, Alliance for Sustainable Energy, LLC, and other contributors. All rights reserved.
+*
+*  Redistribution and use in source and binary forms, with or without modification, are permitted provided that the
+*  following conditions are met:
+*
+*  (1) Redistributions of source code must retain the above copyright notice, this list of conditions and the following
+*  disclaimer.
+*
+*  (2) Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following
+*  disclaimer in the documentation and/or other materials provided with the distribution.
+*
+*  (3) Neither the name of the copyright holder nor the names of any contributors may be used to endorse or promote products
+*  derived from this software without specific prior written permission from the respective party.
+*
+*  (4) Other than as required in clauses (1) and (2), distributions in any form of modifications or other derivative works
+*  may not use the "OpenStudio" trademark, "OS", "os", or any other confusingly similar designation without specific prior
+*  written permission from Alliance for Sustainable Energy, LLC.
+*
+*  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER(S) AND ANY CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
+*  INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+*  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER(S), ANY CONTRIBUTORS, THE UNITED STATES GOVERNMENT, OR THE UNITED
+*  STATES DEPARTMENT OF ENERGY, NOR ANY OF THEIR EMPLOYEES, BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+*  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
+*  USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+*  STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+*  ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+***********************************************************************************************************************/
+
+#include "Polyhedron.hpp"
+#include "Polygon3d.hpp"
+#include "Vector3d.hpp"
+#include "Point3d.hpp"
+#include "Geometry.hpp"
+#include "Plane.hpp"
+#include "Intersection.hpp"
+#include <utilities/geometry/Transformation.hpp>
+
+#include <algorithm>
+#include <utility>
+#include <vector>
+
+namespace openstudio {
+
+Surface3dEdge::Surface3dEdge(Point3d start, Point3d end, Surface3d firstSurface, size_t firstSurfNum)
+  : m_start(std::move(start)), m_end(std::move(end)), m_firstSurfNum(firstSurfNum) {
+  m_allSurfaces.emplace_back(std::move(firstSurface));
+}
+
+/// check equality
+bool Surface3dEdge::operator==(const Surface3dEdge& other) const {
+  return ((isAlmostEqual3dPt(m_start, other.m_start) && isAlmostEqual3dPt(m_end, other.m_end))
+          || (isAlmostEqual3dPt(m_start, other.m_end) && isAlmostEqual3dPt(m_end, other.m_start)));
+}
+
+/// check inequality
+bool Surface3dEdge::operator!=(const Surface3dEdge& other) const {
+  return !(*this == other);
+}
+
+Point3d Surface3dEdge::start() const {
+  return m_start;
+}
+Point3d Surface3dEdge::end() const {
+  return m_end;
+}
+
+size_t Surface3dEdge::count() const {
+  return m_allSurfaces.size();
+}
+
+size_t Surface3dEdge::firstSurfNum() const {
+  return m_firstSurfNum;
+}
+
+std::vector<Surface3d> Surface3dEdge::allSurfaces() const {
+  return m_allSurfaces;
+}
+
+void Surface3dEdge::appendSurface(Surface3d surface) {
+  m_allSurfaces.emplace_back(std::move(surface));
+}
+
+bool Surface3dEdge::containsPoint(const Point3d& testVertex) {
+  return !isAlmostEqual3dPt(m_start, testVertex) && !isAlmostEqual3dPt(m_end, testVertex) && isPointOnLineBetweenPoints(m_start, m_end, testVertex);
+}
+
+std::ostream& operator<<(std::ostream& os, const Surface3dEdge& edge) {
+  os << "Surface3dEdge: start=" << edge.start() << ", end=" << edge.end() << ", count=" << edge.count()
+     << ", firstSurface=" << edge.allSurfaces()[0].name;
+  return os;
+}
+
+Surface3d::Surface3d(std::vector<Point3d> t_vertices, std::string t_name) : vertices(std::move(t_vertices)), name(std::move(t_name)) {}
+
+Polyhedron::Polyhedron(std::vector<Surface3d> surfaces) : m_surfaces(std::move(surfaces)){};
+
+size_t Polyhedron::numVertices() const {
+  size_t count = 0;
+  for (const auto& surface : m_surfaces) {
+    count += surface.vertices.size();
+  }
+  return count;
+}
+
+std::vector<Point3d> Polyhedron::uniqueVertices() const {
+
+  std::vector<Point3d> uniqVertices;
+  uniqVertices.reserve(numVertices());
+
+  for (const auto& surface : m_surfaces) {
+    for (const auto& pt : surface.vertices) {
+      if (std::find_if(uniqVertices.cbegin(), uniqVertices.cend(), [&pt](const auto& unqV) { return isAlmostEqual3dPt(pt, unqV); })
+          == uniqVertices.cend()) {
+        uniqVertices.push_back(pt);
+      }
+    }
+  }
+  return uniqVertices;
+}
+
+std::vector<Surface3dEdge> Polyhedron::edgesNotTwoForEnclosedVolumeTest(const Polyhedron& volumePoly) {
+
+  std::vector<Surface3dEdge> uniqueSurface3dEdges;
+  uniqueSurface3dEdges.reserve(volumePoly.numVertices());
+
+  // construct list of unique edges
+  int surfNum = 0;
+  for (const auto& surface : volumePoly.m_surfaces) {
+    LOG(Debug, "Surface: " << surface.name);
+    const auto& vertices = surface.vertices;
+    for (auto it = vertices.begin(); it != vertices.end(); ++it) {
+      auto itnext = std::next(it);
+      if (itnext == std::end(vertices)) {
+        itnext = std::begin(vertices);
+      }
+      Surface3dEdge thisSurface3dEdge(*it, *itnext, surface, surfNum);
+      auto itFound = std::find(uniqueSurface3dEdges.begin(), uniqueSurface3dEdges.end(), thisSurface3dEdge);
+      if (itFound == uniqueSurface3dEdges.end()) {
+        LOG(Debug, "NOT FOUND: " << thisSurface3dEdge);
+        uniqueSurface3dEdges.emplace_back(std::move(thisSurface3dEdge));
+      } else {
+        LOG(Debug, "    FOUND: " << thisSurface3dEdge);
+        itFound->appendSurface(surface);
+      }
+    }
+    ++surfNum;
+  }
+
+  for (auto& edge : uniqueSurface3dEdges) {
+    LOG(Debug, edge);
+  }
+  // All edges for an enclosed polyhedron should be shared by two (and only two) sides.
+  // So if the count is not two for all edges, the polyhedron is not enclosed, so erase all that are 2
+  uniqueSurface3dEdges.erase(
+    std::remove_if(uniqueSurface3dEdges.begin(), uniqueSurface3dEdges.end(), [](const auto& edge) { return edge.count() == 2; }),
+    uniqueSurface3dEdges.end());
+
+  return uniqueSurface3dEdges;
+}
+
+Polyhedron Polyhedron::updateZonePolygonsForMissingColinearPoints() const {
+  // Make a copy, we don't want to mutate in place
+  Polyhedron updZonePoly(*this);
+
+  std::vector<Point3d> uniqVertices = uniqueVertices();
+
+  for (auto& surface : updZonePoly.m_surfaces) {
+    // for (int iterationLimiter = 0; iterationLimiter < 20; ++iterationLimiter) {  // could probably be while loop but want to make sure it does not get stuck
+    LOG(Debug, surface.name)
+    bool insertedVertext = true;
+    while (insertedVertext) {
+      insertedVertext = false;
+      auto& vertices = surface.vertices;
+      for (auto it = vertices.begin(); it != vertices.end(); ++it) {
+
+        auto itnext = std::next(it);
+        if (itnext == std::end(vertices)) {
+          itnext = std::begin(vertices);
+        }
+
+        // Don't care about surfNum
+        Surface3dEdge thisSurface3dEdge(*it, *itnext, surface, 0);
+
+        // now go through all the vertices and see if they are colinear with start and end vertices
+        for (const auto& testVertex : uniqVertices) {
+          if (thisSurface3dEdge.containsPoint(testVertex)) {
+            LOG(Debug, testVertex << " is on " << thisSurface3dEdge);
+            vertices.insert(itnext, testVertex);
+            insertedVertext = true;
+            break;
+          }
+        }
+        // Break out of the surface loop too, and start again at the while
+        if (insertedVertext) {
+          break;
+        }
+      }
+    }
+  }
+  return updZonePoly;
+}
+
+std::vector<Surface3dEdge> edgesInBoth(const std::vector<Surface3dEdge>& edges1, const std::vector<Surface3dEdge>& edges2) {
+  // this is not optimized but the number of edges for a typical polyhedron is 12 and is probably rarely bigger than 20.
+
+  // (std::set_union could be used, but that requires sorting...)
+  std::vector<Surface3dEdge> inBoth;
+  for (const auto& e1 : edges1) {
+    for (const auto& e2 : edges2) {
+      // If they have the same surface and are the same (within tolerance)
+      if ((e1.firstSurfNum() == e2.firstSurfNum()) && (e1 == e2)) {
+        inBoth.push_back(e1);
+        break;
+      }
+    }
+  }
+  return inBoth;
+}
+
+VolumeEnclosedReturnType Polyhedron::isEnclosedVolume() const {
+
+  VolumeEnclosedReturnType result;
+
+  std::vector<Surface3dEdge> edgeNot2orig = edgesNotTwoForEnclosedVolumeTest(*this);
+  // if all edges had two counts then it is fully enclosed
+  if (edgeNot2orig.empty()) {
+    result.isEnclosedVolume = true;
+    return result;
+  } else {  // if the count is three or greater it is likely that a vertex that is colinear was counted on the faces on one edge and not
+            // on the "other side" of the edge Go through all the points looking for the number that are colinear and see if that is
+            // consistent with the number of edges found that didn't have a count of two
+    Polyhedron updatedZonePoly =
+      updateZonePolygonsForMissingColinearPoints();  // this is done after initial test since it is computationally intensive.
+    std::vector<Surface3dEdge> edgeNot2again = edgesNotTwoForEnclosedVolumeTest(updatedZonePoly);
+    if (edgeNot2again.empty()) {
+      result.isEnclosedVolume = true;
+      return result;
+    }
+    // only return a list of those edges that appear in both the original edge and the revised edges:
+    // this eliminates added edges that will confuse users (edges that were caught by the updateZonePolygonsForMissingColinearPoints routine)
+    result.isEnclosedVolume = false;
+    result.edgesNot2 = edgesInBoth(edgeNot2orig, edgeNot2again);
+    return result;
+  }
+}
+
+// boost::optional<double> Polyhedron::calculatedVolume() const {
+//   auto [isVolEnclosed, edgesNot2] = isEnclosedVolume();
+//   if (isVolEnclosed) {
+//     return calcPolyhedronVolume();
+//   }
+//   return boost::none;
+// }
+
+double Polyhedron::calcPolyhedronVolume() const {
+
+  double volume = 0.0;
+  for (const auto& surface : m_surfaces) {
+    const std::vector<Point3d>& vertices = surface.vertices;
+    Vector3d p3FaceOrigin = vertices[1] - Point3d{};
+    boost::optional<Vector3d> newellAreaVector = getNewellVector(vertices);
+    if (!newellAreaVector) {
+      LOG_AND_THROW("Cannot compute newellAreaVector");
+    }
+    double pyramidVolume = newellAreaVector->dot(p3FaceOrigin);
+    volume += pyramidVolume;
+  }
+  // for (int NFace = 1; NFace <= Poly.NumSurfaceFaces; ++NFace) {
+  //     p3FaceOrigin = Poly.SurfaceFace(NFace).FacePoints(2);
+  //     double pyramidVolume = dot(Poly.SurfaceFace(NFace).NewellAreaVector, (p3FaceOrigin - state.dataVectors->p0));
+  //     Volume += PyramidVolume / 3.0;
+  // }
+  return volume / 6.0;  // Our newellArea vector has twice the length
+}
+
+double Polyhedron::calcDivergenceTheoremVolume() const {
+  double volume = 0.0;
+  for (const auto& surface : m_surfaces) {
+    const std::vector<Point3d>& vertices = surface.vertices;
+    Plane plane(vertices);
+    boost::optional<double> area = getArea(vertices);
+    if (!area) {
+      LOG_AND_THROW("Cannot compute area for " << surface.name);
+    }
+    volume -= plane.d() * area.get() / 3.0;
+  }
+  return volume;
+}
+
+}  // namespace openstudio
diff --git a/src/utilities/geometry/Polyhedron.hpp b/src/utilities/geometry/Polyhedron.hpp
new file mode 100644
index 00000000000..6c12836e2a8
--- /dev/null
+++ b/src/utilities/geometry/Polyhedron.hpp
@@ -0,0 +1,142 @@
+/***********************************************************************************************************************
+*  OpenStudio(R), Copyright (c) 2008-2022, Alliance for Sustainable Energy, LLC, and other contributors. All rights reserved.
+*
+*  Redistribution and use in source and binary forms, with or without modification, are permitted provided that the
+*  following conditions are met:
+*
+*  (1) Redistributions of source code must retain the above copyright notice, this list of conditions and the following
+*  disclaimer.
+*
+*  (2) Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following
+*  disclaimer in the documentation and/or other materials provided with the distribution.
+*
+*  (3) Neither the name of the copyright holder nor the names of any contributors may be used to endorse or promote products
+*  derived from this software without specific prior written permission from the respective party.
+*
+*  (4) Other than as required in clauses (1) and (2), distributions in any form of modifications or other derivative works
+*  may not use the "OpenStudio" trademark, "OS", "os", or any other confusingly similar designation without specific prior
+*  written permission from Alliance for Sustainable Energy, LLC.
+*
+*  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER(S) AND ANY CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
+*  INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+*  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER(S), ANY CONTRIBUTORS, THE UNITED STATES GOVERNMENT, OR THE UNITED
+*  STATES DEPARTMENT OF ENERGY, NOR ANY OF THEIR EMPLOYEES, BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+*  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
+*  USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+*  STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+*  ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+***********************************************************************************************************************/
+
+#ifndef UTILITIES_GEOMETRY_POLYHEDRON_HPP
+#define UTILITIES_GEOMETRY_POLYHEDRON_HPP
+
+#include "../UtilitiesAPI.hpp"
+#include "Point3d.hpp"
+
+#include <vector>
+#include <boost/optional.hpp>
+
+namespace openstudio {
+
+namespace model {
+  class Space;
+  class ThermalZone;
+}  // namespace model
+
+// A thin wrapper to improve reporting, by attaching a name to the vertices
+class UTILITIES_API Surface3d
+{
+ public:
+  Surface3d() = default;
+  Surface3d(std::vector<Point3d> t_vertices, std::string t_name = "None");
+  std::vector<Point3d> vertices;
+  std::string name;
+};
+
+class UTILITIES_API Surface3dEdge
+{
+ public:
+  Surface3dEdge(Point3d start, Point3d end, Surface3d firstSurface, size_t firstSurfNum);
+
+  Point3d start() const;
+  Point3d end() const;
+
+  /// check equality: this uses a tolerance
+  bool operator==(const Surface3dEdge& other) const;
+
+  /// check inequality
+  bool operator!=(const Surface3dEdge& other) const;
+
+  size_t count() const;
+
+  size_t firstSurfNum() const;
+  void appendSurface(Surface3d surface);
+
+  // Checks whether a Point: is not almost equal to the start and end points, and that isPointOnLineBetweenPoints(start, end, testVertex) is true
+  bool containsPoint(const Point3d& testVertex);
+
+  std::vector<Surface3d> allSurfaces() const;
+
+ private:
+  Point3d m_start;
+  Point3d m_end;
+  size_t m_firstSurfNum;
+  std::vector<Surface3d> m_allSurfaces;
+};
+
+using Surface3dEdgeVector = std::vector<Surface3dEdge>;
+
+class VolumeEnclosedReturnType
+{
+ public:
+  VolumeEnclosedReturnType() = default;
+
+  bool isEnclosedVolume = false;
+  Surface3dEdgeVector edgesNot2;
+};
+
+class UTILITIES_API Polyhedron
+{
+ public:
+  Polyhedron() = default;  // Default ctor, for swig
+
+  Polyhedron(std::vector<Surface3d> surfaces);
+
+  // test if the volume described by the polyhedron if full enclosed (would not leak)
+  // This is done by checking that every edge is used exactly TWICE.
+  VolumeEnclosedReturnType isEnclosedVolume() const;
+
+  // Maybe unused
+  std::vector<Point3d> uniqueVertices() const;
+
+  static std::vector<Surface3dEdge> edgesNotTwoForEnclosedVolumeTest(const Polyhedron& volumePoly);
+
+  Polyhedron updateZonePolygonsForMissingColinearPoints() const;
+
+  size_t numVertices() const;
+
+  // Call calculatedVolume as it will handle the logic to dispatch to inidivual methods based on checks
+  // Note JM 2022-05-23: Actually a lot of the needed logic is in PlanarSurface,
+  // so I don't think it's worth trying to fallback to the other calculation modes E+ does (like checking if floor and roof are same and use that etc)
+  // double calculatedVolume() const;
+
+  // Polyhedron MUST be enclosed
+  double calcPolyhedronVolume() const;
+  // Polyhedron MUST be enclosed
+  double calcDivergenceTheoremVolume() const;
+
+ private:
+  REGISTER_LOGGER("utilities.Polyhedron");
+
+  std::vector<Surface3d> m_surfaces;
+};
+
+using PolyhedronVector = std::vector<Polyhedron>;
+using Surface3dVector = std::vector<Surface3d>;
+
+/// ostream operator
+UTILITIES_API std::ostream& operator<<(std::ostream& os, const Surface3dEdge& edge);
+
+}  // namespace openstudio
+
+#endif  // UTILITIES_GEOMETRY_POLYHEDRON_HPP
diff --git a/src/utilities/geometry/Test/Polyhedron_GTest.cpp b/src/utilities/geometry/Test/Polyhedron_GTest.cpp
new file mode 100644
index 00000000000..6000a164dac
--- /dev/null
+++ b/src/utilities/geometry/Test/Polyhedron_GTest.cpp
@@ -0,0 +1,145 @@
+/***********************************************************************************************************************
+*  OpenStudio(R), Copyright (c) 2008-2022, Alliance for Sustainable Energy, LLC, and other contributors. All rights reserved.
+*
+*  Redistribution and use in source and binary forms, with or without modification, are permitted provided that the
+*  following conditions are met:
+*
+*  (1) Redistributions of source code must retain the above copyright notice, this list of conditions and the following
+*  disclaimer.
+*
+*  (2) Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following
+*  disclaimer in the documentation and/or other materials provided with the distribution.
+*
+*  (3) Neither the name of the copyright holder nor the names of any contributors may be used to endorse or promote products
+*  derived from this software without specific prior written permission from the respective party.
+*
+*  (4) Other than as required in clauses (1) and (2), distributions in any form of modifications or other derivative works
+*  may not use the "OpenStudio" trademark, "OS", "os", or any other confusingly similar designation without specific prior
+*  written permission from Alliance for Sustainable Energy, LLC.
+*
+*  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDER(S) AND ANY CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
+*  INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+*  DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER(S), ANY CONTRIBUTORS, THE UNITED STATES GOVERNMENT, OR THE UNITED
+*  STATES DEPARTMENT OF ENERGY, NOR ANY OF THEIR EMPLOYEES, BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+*  EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF
+*  USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
+*  STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
+*  ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+***********************************************************************************************************************/
+
+#include <gtest/gtest.h>
+#include "GeometryFixture.hpp"
+
+#include "../Geometry.hpp"
+#include "../Point3d.hpp"
+#include "../Polyhedron.hpp"
+#include "../PointLatLon.hpp"
+#include "../Vector3d.hpp"
+
+#include <vector>
+
+using namespace openstudio;
+
+std::vector<Point3d> makeRectangle(double min_x = 0.0, double max_x = 10.0, double min_y = 0.0, double max_y = 10.0, double min_z = 0.0,
+                                   double max_z = 0.0) {
+  return {
+    {min_x, min_y, max_z},
+    {min_x, max_y, min_z},
+    {max_x, max_y, min_z},
+    {max_x, min_y, max_z},
+  };
+}
+
+TEST_F(GeometryFixture, Polyhedron_Enclosed) {
+
+  // This is a shoebox model of a lengthy plenum-like zone. 30x10x0.3m. The south wall is split in the middle
+
+  Surface3d south1({{+0.0, +0.0, +0.3}, {+0.0, +0.0, +0.0}, {+15.0, +0.0, +0.0}, {+15.0, +0.0, +0.3}}, "1-SOUTH-1");
+  Surface3d south2({{+15.0, +0.0, +0.3}, {+15.0, +0.0, +0.0}, {+30.0, +0.0, +0.0}, {+30.0, +0.0, +0.3}}, "1-SOUTH-2");
+
+  Surface3d north({{+30.0, +10.0, +0.3}, {+30.0, +10.0, +0.0}, {+0.0, +10.0, +0.0}, {+0.0, +10.0, +0.3}}, "4-NORTH");
+
+  Surface3d east({{+30.0, +0.0, +0.3}, {+30.0, +0.0, +0.0}, {+30.0, +10.0, +0.0}, {+30.0, +10.0, +0.3}}, "3-EAST");
+
+  Surface3d west({{+0.0, +10.0, +0.3}, {+0.0, +10.0, +0.0}, {+0.0, +0.0, +0.0}, {+0.0, +0.0, +0.3}}, "2-WEST");
+
+  Surface3d roof({{+30.0, +0.0, +0.3}, {+30.0, +10.0, +0.3}, {+0.0, +10.0, +0.3}, {+0.0, +0.0, +0.3}}, "ROOF");
+
+  Surface3d floor({{+0.0, +0.0, +0.0}, {+0.0, +10.0, +0.0}, {+30.0, +10.0, +0.0}, {+30.0, +0.0, +0.0}}, "FLOOR");
+
+  Polyhedron zonePoly({south1, south2, north, east, west, roof, floor});
+
+  {
+    // Test individual components
+    // 6 faces for a box, 7 in this case since one face is split in two, each with 4 vertices =>  4*7 = 28
+    EXPECT_EQ(28, zonePoly.numVertices());
+
+    // 8 for the box, plus two for the split of the south wall
+    auto uniqVertices = zonePoly.uniqueVertices();
+    EXPECT_EQ(10, uniqVertices.size());
+
+    std::vector<Surface3dEdge> edgeNot2orig = Polyhedron::edgesNotTwoForEnclosedVolumeTest(zonePoly);
+    EXPECT_EQ(6, edgeNot2orig.size());
+
+    auto updatedZonePoly = zonePoly.updateZonePolygonsForMissingColinearPoints();
+    // We expect the two points from the split south walls to have been added to the roof and floor
+    EXPECT_EQ(30, updatedZonePoly.numVertices());
+
+    std::vector<Surface3dEdge> edgeNot2again = Polyhedron::edgesNotTwoForEnclosedVolumeTest(updatedZonePoly);
+    EXPECT_TRUE(edgeNot2again.empty());
+  }
+
+  // Now do it in one go
+
+  auto r = zonePoly.isEnclosedVolume();
+  EXPECT_TRUE(r.isEnclosedVolume);
+  EXPECT_TRUE(r.edgesNot2.empty());
+
+  double volume = 30.0 * 10.0 * 0.3;
+  EXPECT_EQ(volume, zonePoly.calcPolyhedronVolume());
+  EXPECT_EQ(volume, zonePoly.calcDivergenceTheoremVolume());
+}
+
+TEST_F(GeometryFixture, Polyhedron_Titled_Roof) {
+
+  // This is a 30x10x0.3 base, with a rectangle triangle on top of 30x10x10
+  //                       ▲ z
+  //                       │
+  //                      x├─ 10.0
+  //                    x  │
+  //                  x    │
+  //                x      │
+  //              x        │
+  //            x          │
+  //          x            │
+  //        x              ├─ 0.3
+  //        │              │
+  //   ◄────┼──────────────┼─
+  //  y    10.0            0.0
+
+  Surface3d south2({{+15.0, +0.0, +10.3}, {+15.0, +0.0, +0.0}, {+30.0, +0.0, +0.0}, {+30.0, +0.0, +10.3}}, "1-SOUTH-2");
+
+  // We put extra vertices here to skew the calculate that Space::volume does
+  Surface3d roof({{+30.0, +0.0, +10.3}, {+30.0, +10.0, +0.3}, {+0.0, +10.0, +0.3}, {+0.0, +0.0, +10.3}, {+10.0, +0.0, +10.3}, {+20.0, +0.0, +10.3}},
+                 "ROOF");
+
+  Surface3d east({{+30.0, +0.0, +10.3}, {+30.0, +0.0, +0.0}, {+30.0, +10.0, +0.0}, {+30.0, +10.0, +0.3}}, "3-EAST");
+
+  Surface3d north({{+30.0, +10.0, +0.3}, {+30.0, +10.0, +0.0}, {+0.0, +10.0, +0.0}, {+0.0, +10.0, +0.3}}, "4-NORTH");
+
+  Surface3d west({{+0.0, +10.0, +0.3}, {+0.0, +10.0, +0.0}, {+0.0, +0.0, +0.0}, {+0.0, +0.0, +10.3}}, "2-WEST");
+
+  Surface3d south1({{+0.0, +0.0, +10.3}, {+0.0, +0.0, +0.0}, {+15.0, +0.0, +0.0}, {+15.0, +0.0, +10.3}}, "1-SOUTH-1");
+
+  Surface3d floor({{+0.0, +0.0, +0.0}, {+0.0, +10.0, +0.0}, {+30.0, +10.0, +0.0}, {+30.0, +0.0, +0.0}}, "FLOOR");
+
+  Polyhedron zonePoly({south1, south2, north, east, west, roof, floor});
+
+  auto r = zonePoly.isEnclosedVolume();
+  EXPECT_TRUE(r.isEnclosedVolume);
+  EXPECT_TRUE(r.edgesNot2.empty());
+
+  double volume = 30.0 * 10.0 * 0.3 + 30.0 * 10.0 * 10.0 / 2.0;
+  EXPECT_DOUBLE_EQ(volume, zonePoly.calcPolyhedronVolume());
+  EXPECT_DOUBLE_EQ(volume, zonePoly.calcDivergenceTheoremVolume());
+}


Test patch:
diff --git a/src/model/test/Space_GTest.cpp b/src/model/test/Space_GTest.cpp
index fad604d1907..5b95a242c36 100644
--- a/src/model/test/Space_GTest.cpp
+++ b/src/model/test/Space_GTest.cpp
@@ -3005,4 +3005,60 @@ TEST_F(ModelFixture, DISABLED_ShatteredModel_Existing_3424) {
   EXPECT_EQ(3, ceilingSpace1.size());
   EXPECT_EQ(7, ceilingSpace4.size());
 }
+
+TEST_F(ModelFixture, Space_Polyhedron_Volume) {
+
+  Model m;
+  Space s(m);
+
+  // This is a 30x10x0.3 base, with a rectangle triangle on top of 30x10x10
+  //                       ▲ z
+  //                       │
+  //                      x├─ 10.0
+  //                    x  │
+  //                  x    │
+  //                x      │
+  //              x        │
+  //            x          │
+  //          x            │
+  //        x              ├─ 0.3
+  //        │              │
+  //   ◄────┼──────────────┼─
+  //  y    10.0            0.0
+
+  Surface south2({{+15.0, +0.0, +10.3}, {+15.0, +0.0, +0.0}, {+30.0, +0.0, +0.0}, {+30.0, +0.0, +10.3}}, m);
+  south2.setName("1-SOUTH-2");
+  south2.setSpace(s);
+
+  // Putting extra vertices here on purpose to show that the Space::volume will miscalculate due to averaging foor and ceiling heights
+  Surface roof({{+30.0, +0.0, +10.3}, {+30.0, +10.0, +0.3}, {+0.0, +10.0, +0.3}, {+0.0, +0.0, +10.3}, {+10.0, +0.0, +10.3}, {+20.0, +0.0, +10.3}}, m);
+  roof.setName("ROOF");
+  roof.setSpace(s);
+
+  Surface east({{+30.0, +0.0, +10.3}, {+30.0, +0.0, +0.0}, {+30.0, +10.0, +0.0}, {+30.0, +10.0, +0.3}}, m);
+  east.setName("3-EAST");
+  east.setSpace(s);
+
+  Surface north({{+30.0, +10.0, +0.3}, {+30.0, +10.0, +0.0}, {+0.0, +10.0, +0.0}, {+0.0, +10.0, +0.3}}, m);
+  north.setName("4-NORTH");
+  north.setSpace(s);
+
+  Surface west({{+0.0, +10.0, +0.3}, {+0.0, +10.0, +0.0}, {+0.0, +0.0, +0.0}, {+0.0, +0.0, +10.3}}, m);
+  west.setName("2-WEST");
+  west.setSpace(s);
+
+  Surface south1({{+0.0, +0.0, +10.3}, {+0.0, +0.0, +0.0}, {+15.0, +0.0, +0.0}, {+15.0, +0.0, +10.3}}, m);
+  south1.setName("1-SOUTH-1");
+  south1.setSpace(s);
+
+  Surface floor({{+0.0, +0.0, +0.0}, {+0.0, +10.0, +0.0}, {+30.0, +10.0, +0.0}, {+30.0, +0.0, +0.0}}, m);
+  floor.setName("FLOOR");
+  floor.setSpace(s);
+
+  EXPECT_TRUE(s.isEnclosedVolume());
+
+  double volume = 30.0 * 10.0 * 0.3 + 30.0 * 10.0 * 10.0 / 2.0;
+  EXPECT_EQ(volume, s.volume());
+}
+
 //#  endif // SURFACESHATTERING