Gerard/4372 intersection issue found using create bar

resources/CMakeLists.txt

@@ -119,6 +119,7 @@ set(model_resources_src
   model/floorplan_school.osm
   model/7-7_Windows_Complete.osm
   model/15023_Model12.osm
+  model/offset_tests.osm
 )
 
 

src/model/Surface.cpp

@@ -1458,77 +1458,64 @@ namespace model {
     // 2 - Add all the subsurface rough openings together (so overlap areas don't count twice)
     double Surface_Impl::totalAreaOfSubSurfaces() const {
       double tol = 0.01;
-      // iterate over all sub surfaces
-      // make a map of sub surface/roughOpening vertices (flattened)
-      // iterate over map
-      // check each subsurface for overlap with parent and not overlap with sibling
-      // if either fails, replace roughOpening vertices with original vertices
-      // Add up all the openign areas
-
-      // Records the rough opening of each sub surface
-      std::map<SubSurface, Point3dVector> roughOpenings;
-      // Get the flattened parent surface vertices
-      Transformation parentToXY = Transformation::alignFace(this->vertices()).inverse();
-      std::vector<Point3d> parentVertices = parentToXY * this->vertices();
-      // Make sure the parent surface is oriented clockwise
-      auto norm = openstudio::getOutwardNormal(parentVertices);
+
+      // There must be a simpler way - after all we are only intersetced in areas
+      // 1 - Get the surface polygon flattened
+      // 2 - Get all the sub-surface polygons with frame offset applied
+      // 3 - STart with first poly
+      // 4 - Does it overlap any other polys? (A & B)
+      // 5 - Y: Get the combined area of the two polygons, remove polygon A & B from the list insert polygon C
+      //     N: Remove polygon A form the list add to finished polygon list
+      // 6 When no more overlaps are detected we are finished. Get the total area of all the polygons
+
+      // Simple case - no sub-surfaces
+      if (this->subSurfaces().size() == 0) return 0;
+
+      // Flattened surface vertices
+      Transformation surfaceToXY = Transformation::alignFace(this->vertices()).inverse();
+      std::vector<Point3d> surfaceVertices = surfaceToXY * this->vertices();
+      // Make sure the surface surface is oriented clockwise
+      auto norm = openstudio::getOutwardNormal(surfaceVertices);
       if (norm && norm->z() > 0) {
-        std::reverse(parentVertices.begin(), parentVertices.end());
+        std::reverse(surfaceVertices.begin(), surfaceVertices.end());
       }
 
+      Point3dVectorVector subSurfaces;
+
       // Get the flattened sub surface vertices for windows
       for (const SubSurface& subSurface : this->subSurfaces()) {
         if (istringEqual(subSurface.subSurfaceType(), "FixedWindow") || istringEqual(subSurface.subSurfaceType(), "OperableWindow")) {
-          auto opening = parentToXY * subSurface.roughOpeningVertices();
-          auto norm = openstudio::getOutwardNormal(opening);
-          if (norm && norm->z() > 0) {
-            std::reverse(opening.begin(), opening.end());
-          }
-          roughOpenings[subSurface] = opening;
-        }
-      }
-
-      // Check each rough opening against the parent surface, if any vertices are outside
-      // the parent surface then revert back to the original vertices
-      for (const auto& opening : roughOpenings) {
-        const SubSurface& subSurface = opening.first;
-        std::vector<Point3d> flattenedVertices = opening.second;
-
-        if (!openstudio::polygonInPolygon(flattenedVertices, parentVertices, tol)) {
-          // One or more vertices not within the parent surface so replace with the original vertices
-          auto vertices = parentToXY * subSurface.vertices();
-          auto norm = openstudio::getOutwardNormal(vertices);
+          auto roughOpening = surfaceToXY * subSurface.roughOpeningVertices();
+          auto norm = openstudio::getOutwardNormal(roughOpening);
           if (norm && norm->z() > 0) {
-            std::reverse(vertices.begin(), vertices.end());
+            std::reverse(roughOpening.begin(), roughOpening.end());
           }
-          roughOpenings[subSurface] = vertices;
+          subSurfaces.push_back(roughOpening);
         }
       }
 
-      // Check subsurfaces for overlaps
-      // (if I knew how I'd order themn left to right in x then we could test adjacent ones only)
-      for (const auto& opening1 : roughOpenings) {
-        for (const auto& opening2 : roughOpenings) {
-          if (opening1 == opening2) {
-            continue;
-          }
+      // Add all the subsurfaces togther
+      auto finalSurfaces = openstudio::joinAll(subSurfaces, tol);
 
-          if (auto result = openstudio::intersect(opening1.second, opening2.second, tol)) {
-            //There's an overlap so swap out the overlapped vertices for the original
-            auto vertices = parentToXY * opening1.first.vertices();
-            auto norm = openstudio::getOutwardNormal(vertices);
-            if (norm && norm->z() > 0) {
-              std::reverse(vertices.begin(), vertices.end());
+      double area = 0;
+      // No overlappingsurfaces so we can add up the surface areas including multipliers
+      if (finalSurfaces.size() == subSurfaces.size()) {
+        for (const auto& subSurface : this->subSurfaces()) {
+          if (istringEqual(subSurface.subSurfaceType(), "FixedWindow") || istringEqual(subSurface.subSurfaceType(), "OperableWindow")) {
+            auto roughOpening = surfaceToXY * subSurface.roughOpeningVertices();
+            if (!openstudio::polygonInPolygon(roughOpening, surfaceVertices, tol)) {
+              roughOpening = surfaceToXY * subSurface.vertices();
             }
-            roughOpenings[opening1.first] = vertices;
+            area += openstudio::getArea(roughOpening).value() * subSurface.multiplier();
           }
         }
+
+        return area;
       }
 
-      // Accumulate the areas
-      double area = std::accumulate(roughOpenings.cbegin(), roughOpenings.cend(), 0.0, [](double sum, auto& roughOpening) {
-        return sum + openstudio::getArea(roughOpening.second).value() * roughOpening.first.multiplier();
-      });
+      for (const auto& subSurface : subSurfaces) {
+        area += openstudio::getArea(subSurface).value();
+      }
 
       return area;
     }

src/model/test/Model_GTest.cpp

@@ -82,6 +82,8 @@
 #include "../../utilities/idf/WorkspaceObject.hpp"
 #include "../../utilities/idf/ValidityReport.hpp"
 
+#include "../../osversion/VersionTranslator.hpp"
+
 #include <utilities/idd/IddEnums.hxx>
 
 #include <boost/algorithm/string/case_conv.hpp>
@@ -790,3 +792,59 @@ TEST_F(ModelFixture, Ensure_Name_Unicity_SpaceAndSpaceGroupNames) {
     EXPECT_NE(mos[i1].nameString(), mos[i2].nameString());
   }
 }
+
+TEST_F(ModelFixture, Issue_4372) {
+
+  // Open OffsetTests
+  // INtersect surfaces
+  // Match Surfaces
+  // Check the internal walls are paired
+  // Surface 4 Space 101 with SPace 102 surface 12
+  // Surface 2 Space 102 with
+
+  openstudio::path modelPath = resourcesPath() / "model" / toPath("offset_tests.osm");
+  ASSERT_TRUE(openstudio::filesystem::exists(modelPath));
+
+  openstudio::osversion::VersionTranslator vt;
+  boost::optional<openstudio::model::Model> model = vt.loadModel(modelPath);
+  ASSERT_TRUE(model);
+
+  std::vector<Space> spaces = model->getConcreteModelObjects<Space>();
+  intersectSurfaces(spaces);
+
+  matchSurfaces(spaces);
+  std::vector<Surface> surfacesAfter = model->getConcreteModelObjects<Surface>();
+  for (const auto& surface : surfacesAfter) {
+
+    std::string name = surface.name().value();
+    OptionalSurface otherSurface;
+    if (name == "Surface 4") {
+      otherSurface = surface.adjacentSurface();
+      ASSERT_TRUE(otherSurface);
+      EXPECT_EQ(otherSurface->nameString(), "Surface 8");
+    } else if (name == "Surface 10") {
+      otherSurface = surface.adjacentSurface();
+      ASSERT_TRUE(otherSurface);
+      ASSERT_EQ(otherSurface->name().value(), "Surface 14");
+    } else if (name == "Surface 16") {
+      otherSurface = surface.adjacentSurface();
+      ASSERT_TRUE(otherSurface);
+      ASSERT_EQ(otherSurface->name().value(), "Surface 20");
+    } else if (name == "Surface 22") {
+      otherSurface = surface.adjacentSurface();
+      ASSERT_TRUE(otherSurface);
+      ASSERT_EQ(otherSurface->name().value(), "Surface 26");
+    } else if (name == "Surface 28") {
+      otherSurface = surface.adjacentSurface();
+      ASSERT_TRUE(otherSurface);
+      ASSERT_EQ(otherSurface->name().value(), "Surface 32");
+    }
+
+    // if (otherSurface) {
+    //   LOG(Info, "Surface " << surface.name().value() << " is paired with " << otherSurface->name().value());
+    // }
+  }
+
+  // modelPath = resourcesPath() / "model" / toPath("offset_tests_matched.osm");
+  // model->save(modelPath, true);
+}

src/model/test/Space_GTest.cpp

@@ -1716,6 +1716,7 @@ TEST_F(ModelFixture, Space_intersectSurfaces_degenerate2) {
   std::vector<Point3d> vertices;
 
   // bottom floor
+  BuildingStory bottomStory(m);
 
   // bottom core
   vertices.clear();
@@ -1726,6 +1727,7 @@ TEST_F(ModelFixture, Space_intersectSurfaces_degenerate2) {
   boost::optional<Space> bottomCore = Space::fromFloorPrint(vertices, 3, m);
   ASSERT_TRUE(bottomCore);
   bottomCore->setZOrigin(0);
+  bottomCore->setBuildingStory(bottomStory);
 
   // bottom top
   vertices.clear();
@@ -1736,6 +1738,7 @@ TEST_F(ModelFixture, Space_intersectSurfaces_degenerate2) {
   boost::optional<Space> bottomTop = Space::fromFloorPrint(vertices, 3, m);
   ASSERT_TRUE(bottomTop);
   bottomTop->setZOrigin(0);
+  bottomTop->setBuildingStory(bottomStory);
 
   // bottom right
   vertices.clear();
@@ -1746,6 +1749,7 @@ TEST_F(ModelFixture, Space_intersectSurfaces_degenerate2) {
   boost::optional<Space> bottomRight = Space::fromFloorPrint(vertices, 3, m);
   ASSERT_TRUE(bottomRight);
   bottomRight->setZOrigin(0);
+  bottomRight->setBuildingStory(bottomStory);
 
   // bottom bottom
   vertices.clear();
@@ -1756,6 +1760,7 @@ TEST_F(ModelFixture, Space_intersectSurfaces_degenerate2) {
   boost::optional<Space> bottomBottom = Space::fromFloorPrint(vertices, 3, m);
   ASSERT_TRUE(bottomBottom);
   bottomBottom->setZOrigin(0);
+  bottomBottom->setBuildingStory(bottomStory);
 
   // bottom left
   vertices.clear();
@@ -1766,8 +1771,10 @@ TEST_F(ModelFixture, Space_intersectSurfaces_degenerate2) {
   boost::optional<Space> bottomLeft = Space::fromFloorPrint(vertices, 3, m);
   ASSERT_TRUE(bottomLeft);
   bottomLeft->setZOrigin(0);
+  bottomLeft->setBuildingStory(bottomStory);
 
   // top floor
+  BuildingStory topStory(m);
 
   // top core
   vertices.clear();
@@ -1778,6 +1785,7 @@ TEST_F(ModelFixture, Space_intersectSurfaces_degenerate2) {
   boost::optional<Space> topCore = Space::fromFloorPrint(vertices, 3, m);
   ASSERT_TRUE(topCore);
   topCore->setZOrigin(3);
+  topCore->setBuildingStory(topStory);
 
   // top top
   vertices.clear();
@@ -1788,6 +1796,7 @@ TEST_F(ModelFixture, Space_intersectSurfaces_degenerate2) {
   boost::optional<Space> topTop = Space::fromFloorPrint(vertices, 3, m);
   ASSERT_TRUE(topTop);
   topTop->setZOrigin(3);
+  topTop->setBuildingStory(topStory);
 
   // top right
   vertices.clear();
@@ -1798,6 +1807,7 @@ TEST_F(ModelFixture, Space_intersectSurfaces_degenerate2) {
   boost::optional<Space> topRight = Space::fromFloorPrint(vertices, 3, m);
   ASSERT_TRUE(topRight);
   topRight->setZOrigin(3);
+  topRight->setBuildingStory(topStory);
 
   // top bottom
   vertices.clear();
@@ -1808,6 +1818,7 @@ TEST_F(ModelFixture, Space_intersectSurfaces_degenerate2) {
   boost::optional<Space> topBottom = Space::fromFloorPrint(vertices, 3, m);
   ASSERT_TRUE(topBottom);
   topBottom->setZOrigin(3);
+  topBottom->setBuildingStory(topStory);
 
   // top left
   vertices.clear();
@@ -1818,14 +1829,20 @@ TEST_F(ModelFixture, Space_intersectSurfaces_degenerate2) {
   boost::optional<Space> topLeft = Space::fromFloorPrint(vertices, 3, m);
   ASSERT_TRUE(topLeft);
   topLeft->setZOrigin(3);
+  topLeft->setBuildingStory(topStory);
 
   // create thermal zones
   std::vector<Space> spaces = m.getConcreteModelObjects<Space>();
+  //std::sort(spaces.begin(), spaces.end(), WorkspaceObjectNameLess());
+
   for (auto& space : spaces) {
     ThermalZone z(m);
     space.setThermalZone(z);
   }
 
+  //openstudio::path outpath1 = resourcesPath() / toPath("model/Space_intersectSurfaces_degenerate2_before_intersect.osm");
+  //m.save(outpath1, true);
+
   intersectSurfaces(spaces);
   matchSurfaces(spaces);
 
@@ -1838,6 +1855,8 @@ TEST_F(ModelFixture, Space_intersectSurfaces_degenerate2) {
   int numRoofSurfaces = 0;
 
   std::vector<Surface> surfaces = m.getConcreteModelObjects<Surface>();
+  //std::sort(surfaces.begin(), surfaces.end(), WorkspaceObjectNameLess());
+
   for (auto& surface : surfaces) {
     if (istringEqual(surface.surfaceType(), "RoofCeiling")) {
       if (istringEqual(surface.outsideBoundaryCondition(), "Outdoors")) {

src/model/test/SubSurface_GTest.cpp

@@ -1453,9 +1453,11 @@ TEST_F(ModelFixture, Issue_4361_Multi_Subsurfaces_Non_Overlapping) {
   subSurface2.setSurface(surface);
   subSurface2.setSubSurfaceType("FixedWindow");
 
+  // Surface area is 8, sub-surface areas are 2 therefore we should expect a WWR of 0.25
   double windowWallRatio = surface.windowToWallRatio();
   EXPECT_NEAR(windowWallRatio, 0.25, 0.01);
 
+  // Set frame and divider on both sub-surfaces
   WindowPropertyFrameAndDivider frame1(model);
   frame1.setFrameWidth(0.030);
   subSurface1.setWindowPropertyFrameAndDivider(frame1);
@@ -1464,6 +1466,7 @@ TEST_F(ModelFixture, Issue_4361_Multi_Subsurfaces_Non_Overlapping) {
   frame2.setFrameWidth(0.030);
   subSurface2.setWindowPropertyFrameAndDivider(frame2);
 
+  // Surface area is still 8, sub-surface areas are 1.06x1.06x2 = 2.2472 so we should expect a WWR of 0.2809 (2.2472 / 8)
   windowWallRatio = surface.windowToWallRatio();
   EXPECT_NEAR(windowWallRatio, 0.281, 0.01);
 }
@@ -1510,9 +1513,11 @@ TEST_F(ModelFixture, Issue_4361_Multi_Subsurfaces_Overlapping) {
   subSurface2.setSurface(surface);
   subSurface2.setSubSurfaceType("FixedWindow");
 
+  // Surface area is 8, sub surface areas are 2 therefore we should expect a WWR of 0.25
   double windowWallRatio = surface.windowToWallRatio();
   EXPECT_NEAR(windowWallRatio, 0.25, 0.01);
 
+  // Set a frame and divider on both sub-surfaces
   WindowPropertyFrameAndDivider frame1(model);
   frame1.setFrameWidth(0.030);
   subSurface1.setWindowPropertyFrameAndDivider(frame1);
@@ -1521,6 +1526,8 @@ TEST_F(ModelFixture, Issue_4361_Multi_Subsurfaces_Overlapping) {
   frame2.setFrameWidth(0.030);
   subSurface2.setWindowPropertyFrameAndDivider(frame2);
 
+  // When the frame and dividers are added the surfaces will overlap to form a single
+  // sub-surface which is 2.07 x 1.06 = 2.1942 which gives a WWR of 0.274 (2.1942 / 8)
   windowWallRatio = surface.windowToWallRatio();
-  EXPECT_NEAR(windowWallRatio, 0.2654, 0.01);
+  EXPECT_NEAR(windowWallRatio, 0.2742, 0.01);
 }

src/utilities/geometry/Intersection.cpp

@@ -70,10 +70,11 @@ std::vector<BoostPolygon> removeSpikesEx(const BoostPolygon& polygon) {
 
   //const double buffer_distance = 1.0;
   //const int points_per_circle = 36;
-  const double amount = 0.005;
+  const double offsetBy = 0.005;
+  const double tol = offsetBy;
 
-  boost::geometry::strategy::buffer::distance_symmetric<coordinate_type> expand(amount);
-  boost::geometry::strategy::buffer::distance_symmetric<coordinate_type> shrink(-amount);
+  boost::geometry::strategy::buffer::distance_symmetric<coordinate_type> expand(offsetBy);
+  boost::geometry::strategy::buffer::distance_symmetric<coordinate_type> shrink(-offsetBy);
   boost::geometry::strategy::buffer::join_miter join_strategy;
   boost::geometry::strategy::buffer::end_flat end_strategy;
   boost::geometry::strategy::buffer::side_straight side_strategy;
@@ -87,11 +88,16 @@ std::vector<BoostPolygon> removeSpikesEx(const BoostPolygon& polygon) {
   BoostMultiPolygon result;
   boost::geometry::buffer(polygon, resultShrink, shrink, side_strategy, join_strategy, end_strategy, point_strategy);
   boost::geometry::buffer(resultShrink, resultExpand, expand, side_strategy, join_strategy, end_strategy, point_strategy);
-  boost::geometry::simplify(resultExpand, result, amount);
+  boost::geometry::simplify(resultExpand, result, offsetBy);
 
   // cppcheck-suppress constStatement
   std::vector<BoostPolygon> solution;
   if (result.size() == 0) {
+    // There's no result so return the original polygon
+    solution.push_back(polygon);
+    return solution;
+  } else if (result.size() == 1 && result[0].outer().size() == polygon.outer().size()) {
+    // Number of vertices didn't change so no spikes were removed
     solution.push_back(polygon);
     return solution;
   } else {
@@ -105,7 +111,7 @@ std::vector<BoostPolygon> removeSpikesEx(const BoostPolygon& polygon) {
         for (unsigned j = 0; j < polygon.outer().size(); j++) {
           Point3d p1(polygon.outer()[j].x(), polygon.outer()[j].y(), 0);
           // Two points are within tolerance set the result to the original input point
-          if (getDistance(point3d, p1) <= 0.05) {
+          if (getDistance(point3d, p1) <= tol) {
             boostPolygon.outer()[i].x(polygon.outer()[j].x());
             boostPolygon.outer()[i].y(polygon.outer()[j].y());
             break;