Fix#127 - Fix volume calculation to match EnergyPlus

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>
 

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

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

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");
 

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

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

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

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

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 {