#3644 Move hex overlap code to RigCellGeometryTools and clamp in z-direction as well.

* The latter is not important for 2d Maps because the height of the 2d extrusion is always at least as much as the cells.
* However this makes it so the code can estimate overlap with bounding boxes that do not reach the full height of the hexahedron.
* Also added unit tests.

ApplicationCode/ProjectDataModel/Rim2dGridProjection.cpp

@@ -684,7 +684,7 @@ std::vector<std::pair<size_t, double>> Rim2dGridProjection::visibleCellsAndWeigh
             localGrid->cellCornerVertices(localCellIdx, hexCorners.data());
 
             cvf::BoundingBox overlapBBox;
-            std::array<cvf::Vec3d, 8> overlapCorners = createHexOverlapEstimation(bbox2dElement, hexCorners, &overlapBBox);
+            std::array<cvf::Vec3d, 8> overlapCorners = RigCellGeometryTools::estimateHexOverlapWithBoundingBox(hexCorners, bbox2dElement, &overlapBBox);
             
             double overlapVolume = RigCellGeometryTools::calculateCellVolume(overlapCorners);
 
@@ -815,41 +815,6 @@ double Rim2dGridProjection::findSoilResult(size_t cellGlobalIdx) const
     return 0.0;
 }
 
-//--------------------------------------------------------------------------------------------------
-///
-//--------------------------------------------------------------------------------------------------
-std::array<cvf::Vec3d, 8> Rim2dGridProjection::createHexOverlapEstimation(const cvf::BoundingBox& bbox2dElement, const std::array<cvf::Vec3d, 8>& hexCorners, cvf::BoundingBox* overlapBoundingBox) const
-{
-    CVF_ASSERT(overlapBoundingBox);
-    *overlapBoundingBox = cvf::BoundingBox();
-    std::array<cvf::Vec3d, 8> overlapCorners = hexCorners;
-    // A reasonable approximation to the overlap volume
-    cvf::Plane topPlane;    topPlane.setFromPoints(hexCorners[0], hexCorners[1], hexCorners[2]);
-    cvf::Plane bottomPlane; bottomPlane.setFromPoints(hexCorners[4], hexCorners[5], hexCorners[6]);
-
-    for (size_t i = 0; i < 4; ++i)
-    {
-        cvf::Vec3d& corner = overlapCorners[i];
-        corner.x() = cvf::Math::clamp(corner.x(), bbox2dElement.min().x(), bbox2dElement.max().x());
-        corner.y() = cvf::Math::clamp(corner.y(), bbox2dElement.min().y(), bbox2dElement.max().y());
-        cvf::Vec3d maxZCorner = corner; maxZCorner.z() = bbox2dElement.max().z();
-        cvf::Vec3d minZCorner = corner; minZCorner.z() = bbox2dElement.min().z();
-        topPlane.intersect(maxZCorner, minZCorner, &corner);
-        overlapBoundingBox->add(corner);
-    }
-    for (size_t i = 4; i < 8; ++i)
-    {
-        cvf::Vec3d& corner = overlapCorners[i];
-        corner.x() = cvf::Math::clamp(corner.x(), bbox2dElement.min().x(), bbox2dElement.max().x());
-        corner.y() = cvf::Math::clamp(corner.y(), bbox2dElement.min().y(), bbox2dElement.max().y());
-        cvf::Vec3d maxZCorner = corner; maxZCorner.z() = bbox2dElement.max().z();
-        cvf::Vec3d minZCorner = corner; minZCorner.z() = bbox2dElement.min().z();
-        bottomPlane.intersect(maxZCorner, minZCorner, &corner);
-        overlapBoundingBox->add(corner);
-    }
-    return overlapCorners;
-}
-
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp

@@ -27,6 +27,8 @@
 
 #include "clipper/clipper.hpp"
 
+#include "cvfMath.h"
+
 #include <vector>
 #include <array>
 
@@ -82,6 +84,43 @@ double RigCellGeometryTools::calculateCellVolume(const std::array<cvf::Vec3d, 8>
     return std::abs(volume); // Altogether 18 + 3*17 + 3 + 1 flops = 73 flops.
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::array<cvf::Vec3d, 8> RigCellGeometryTools::estimateHexOverlapWithBoundingBox(const std::array<cvf::Vec3d, 8>& hexCorners, const cvf::BoundingBox& boundingBox, cvf::BoundingBox* overlapBoundingBox)
+{
+    CVF_ASSERT(overlapBoundingBox);
+    *overlapBoundingBox = cvf::BoundingBox();
+    std::array<cvf::Vec3d, 8> overlapCorners = hexCorners;
+    // A reasonable approximation to the overlap volume
+    cvf::Plane topPlane;    topPlane.setFromPoints(hexCorners[0], hexCorners[1], hexCorners[2]);
+    cvf::Plane bottomPlane; bottomPlane.setFromPoints(hexCorners[4], hexCorners[5], hexCorners[6]);
+
+    for (size_t i = 0; i < 4; ++i)
+    {
+        cvf::Vec3d& corner = overlapCorners[i];
+        corner.x() = cvf::Math::clamp(corner.x(), boundingBox.min().x(), boundingBox.max().x());
+        corner.y() = cvf::Math::clamp(corner.y(), boundingBox.min().y(), boundingBox.max().y());
+        corner.z() = cvf::Math::clamp(corner.z(), boundingBox.min().z(), boundingBox.max().z());
+        cvf::Vec3d maxZCorner = corner; maxZCorner.z() = boundingBox.max().z();
+        cvf::Vec3d minZCorner = corner; minZCorner.z() = boundingBox.min().z();
+        topPlane.intersect(minZCorner, maxZCorner, &corner);
+        overlapBoundingBox->add(corner);
+    }
+    for (size_t i = 4; i < 8; ++i)
+    {
+        cvf::Vec3d& corner = overlapCorners[i];
+        corner.x() = cvf::Math::clamp(corner.x(), boundingBox.min().x(), boundingBox.max().x());
+        corner.y() = cvf::Math::clamp(corner.y(), boundingBox.min().y(), boundingBox.max().y());
+        corner.z() = cvf::Math::clamp(corner.z(), boundingBox.min().z(), boundingBox.max().z());
+        cvf::Vec3d maxZCorner = corner; maxZCorner.z() = boundingBox.max().z();
+        cvf::Vec3d minZCorner = corner; minZCorner.z() = boundingBox.min().z();
+        bottomPlane.intersect(minZCorner, maxZCorner, &corner);
+        overlapBoundingBox->add(corner);
+    }
+    return overlapCorners;
+}
+
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h

@@ -28,16 +28,21 @@
 #include <list>
 #include <vector>
 
-
-
 class RigCellGeometryTools
 {
 public:
     static double                    calculateCellVolume(const std::array<cvf::Vec3d, 8>& hexCorners);
+    static std::array<cvf::Vec3d, 8> estimateHexOverlapWithBoundingBox(const std::array<cvf::Vec3d, 8>& hexCorners,
+                                                                       const cvf::BoundingBox&          boundingBox2dExtrusion,
+                                                                       cvf::BoundingBox*                overlapBoundingBox);
 
-    static void createPolygonFromLineSegments(std::list<std::pair<cvf::Vec3d, cvf::Vec3d>> &intersectionLineSegments, std::vector<std::vector<cvf::Vec3d>> &polygons);
+    static void createPolygonFromLineSegments(std::list<std::pair<cvf::Vec3d, cvf::Vec3d>>& intersectionLineSegments,
+                                              std::vector<std::vector<cvf::Vec3d>>&         polygons);
 
-    static void findCellLocalXYZ(const std::array<cvf::Vec3d, 8>& hexCorners, cvf::Vec3d &localXdirection, cvf::Vec3d &localYdirection, cvf::Vec3d &localZdirection);
+    static void findCellLocalXYZ(const std::array<cvf::Vec3d, 8>& hexCorners,
+                                 cvf::Vec3d&                      localXdirection,
+                                 cvf::Vec3d&                      localYdirection,
+                                 cvf::Vec3d&                      localZdirection);
 
     static double polygonLengthInLocalXdirWeightedByArea(const std::vector<cvf::Vec3d>& polygon2d);
 
@@ -47,19 +52,26 @@ public:
     static std::vector<std::vector<cvf::Vec3d>> subtractPolygons(const std::vector<cvf::Vec3d>&              sourcePolygon,
                                                                  const std::vector<std::vector<cvf::Vec3d>>& polygonToSubtract);
 
-    enum ZInterpolationType { INTERPOLATE_LINE_Z, USE_HUGEVAL, USE_ZERO};
+    enum ZInterpolationType
+    {
+        INTERPOLATE_LINE_Z,
+        USE_HUGEVAL,
+        USE_ZERO
+    };
     static std::vector<std::vector<cvf::Vec3d>> clipPolylineByPolygon(const std::vector<cvf::Vec3d>& polyLine,
                                                                       const std::vector<cvf::Vec3d>& polygon,
                                                                       ZInterpolationType             interpolType = USE_ZERO);
 
-    static std::pair<cvf::Vec3d, cvf::Vec3d> getLineThroughBoundingBox(const cvf::Vec3d& lineDirection, const cvf::BoundingBox& polygonBBox, const cvf::Vec3d& pointOnLine);
+    static std::pair<cvf::Vec3d, cvf::Vec3d> getLineThroughBoundingBox(const cvf::Vec3d&       lineDirection,
+                                                                       const cvf::BoundingBox& polygonBBox,
+                                                                       const cvf::Vec3d&       pointOnLine);
 
-    static double getLengthOfPolygonAlongLine(const std::pair<cvf::Vec3d, cvf::Vec3d>& line, const std::vector<cvf::Vec3d>& polygon);
+    static double getLengthOfPolygonAlongLine(const std::pair<cvf::Vec3d, cvf::Vec3d>& line,
+                                              const std::vector<cvf::Vec3d>&           polygon);
 
     static std::vector<cvf::Vec3d> unionOfPolygons(const std::vector<std::vector<cvf::Vec3d>>& polygons);
 
 private:
     static std::vector<cvf::Vec3d> ajustPolygonToAvoidIntersectionsAtVertex(const std::vector<cvf::Vec3d>& polyLine,
                                                                             const std::vector<cvf::Vec3d>& polygon);
-
 };

ApplicationCode/UnitTests/RigCellGeometryTools-Test.cpp

@@ -26,21 +26,101 @@
 //--------------------------------------------------------------------------------------------------
 TEST(RigCellGeometryTools, calculateCellVolumeTest)
 {
-    cvf::BoundingBox bbox(cvf::Vec3d(1.0, -2.0, 5.0), cvf::Vec3d(200.0, 3.0, 1500.0));
+    cvf::BoundingBox bbox(cvf::Vec3d(1.0, -2.0, 5.0), cvf::Vec3d(500.0, 3.0, 1500.0));
     cvf::Vec3d extent = bbox.extent();
     double bboxVolume = extent.x() * extent.y() * extent.z();
 
     std::array<cvf::Vec3d, 8> cornerVertices;
     bbox.cornerVertices(cornerVertices.data());    
 
+    // This is a cuboid. The result should be exact
     EXPECT_DOUBLE_EQ(bboxVolume, RigCellGeometryTools::calculateCellVolume(cornerVertices));
 
-    cornerVertices[1].x() += 100.0;
-    cornerVertices[2].x() += 100.0;
+    // Distort it by adding a tetrahedron to the volume
+    cornerVertices[1].x() += bbox.extent().x() / 3.0;
+    cornerVertices[2].x() += bbox.extent().x() / 3.0;
 
-    double extraVolume = 0.5 * extent.z() * 100 * extent.y();
+    double extraVolume = 0.5 * extent.z() * bbox.extent().x() / 3.0 * extent.y();
     
     EXPECT_DOUBLE_EQ(bboxVolume + extraVolume, RigCellGeometryTools::calculateCellVolume(cornerVertices));
+
+    // The overlap with the original bounding box should just yield the original bounding box
+    cvf::BoundingBox overlapBoundingBox;
+    std::array<cvf::Vec3d, 8> overlapVertices = RigCellGeometryTools::estimateHexOverlapWithBoundingBox(cornerVertices, bbox, &overlapBoundingBox);
+
+    EXPECT_DOUBLE_EQ(bboxVolume, RigCellGeometryTools::calculateCellVolume(overlapVertices));
+
+    cvf::Vec3d overlapExtent = overlapBoundingBox.extent();
+    double overlapBBoxVolume = overlapExtent.x() * overlapExtent.y() * overlapExtent.z();
+    EXPECT_DOUBLE_EQ(bboxVolume, overlapBBoxVolume);
+    
+    // Shift bounding box by half the original extent in x-direction.
+    // It should now contain the full tetrahedron + half the original bounding box
+    std::array<cvf::Vec3d, 8> tetrahedronBoxVertices;
+    bbox.cornerVertices(tetrahedronBoxVertices.data());
+    cvf::BoundingBox tetrahedronBBox;
+    for (cvf::Vec3d& corner : tetrahedronBoxVertices)
+    {
+        corner.x() += 0.5 * bbox.extent().x();
+        tetrahedronBBox.add(corner);
+    }
+    overlapVertices = RigCellGeometryTools::estimateHexOverlapWithBoundingBox(cornerVertices, tetrahedronBBox, &overlapBoundingBox);
+
+    EXPECT_DOUBLE_EQ(bboxVolume * 0.5 + extraVolume, RigCellGeometryTools::calculateCellVolume(overlapVertices));
+
+    // Shift it the rest of the original extent in x-direction.
+    // The bounding box should now contain only the tetrahedron.
+    tetrahedronBBox = cvf::BoundingBox();
+    for (cvf::Vec3d& corner : tetrahedronBoxVertices)
+    {
+        corner.x() += 0.5 * bbox.extent().x();
+        tetrahedronBBox.add(corner);
+    }
+    overlapVertices = RigCellGeometryTools::estimateHexOverlapWithBoundingBox(cornerVertices, tetrahedronBBox, &overlapBoundingBox);
+
+    EXPECT_DOUBLE_EQ(extraVolume, RigCellGeometryTools::calculateCellVolume(overlapVertices));
+
+    // Expand original bounding box to be much larger than the hex
+    bbox.expand(2000);
+    overlapVertices = RigCellGeometryTools::estimateHexOverlapWithBoundingBox(cornerVertices, bbox, &overlapBoundingBox);
+    EXPECT_DOUBLE_EQ(bboxVolume + extraVolume, RigCellGeometryTools::calculateCellVolume(overlapVertices));    
+}
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+TEST(RigCellGeometryTools, calculateCellVolumeTest2)
+{
+    cvf::BoundingBox bbox(cvf::Vec3d(0.0, 0.0, 0.0), cvf::Vec3d(100.0, 100.0, 100.0));
+    std::array<cvf::Vec3d, 8> cornerVertices;
+    bbox.cornerVertices(cornerVertices.data());
+
+    cornerVertices[5].z() = cornerVertices[1].z();
+    cornerVertices[6].z() = cornerVertices[2].z();
+
+    double totalCellVolume = 0.5 * 100.0 * 100.0 * 100.0;
+    EXPECT_DOUBLE_EQ(totalCellVolume, RigCellGeometryTools::calculateCellVolume(cornerVertices));
+    cvf::BoundingBox innerBBox(cvf::Vec3d(25.0, 25.0, -10.0), cvf::Vec3d(75.0, 75.0, 110.0));
+
+    double expectedOverlap = 50 * 50 * 25 + 0.5 * 50 * 50 * 50;
+
+    cvf::BoundingBox overlapBoundingBox;
+    std::array<cvf::Vec3d, 8> overlapVertices = RigCellGeometryTools::estimateHexOverlapWithBoundingBox(cornerVertices, innerBBox, &overlapBoundingBox);
+    EXPECT_DOUBLE_EQ(expectedOverlap, RigCellGeometryTools::calculateCellVolume(overlapVertices));
+
+    cvf::BoundingBox smallerInnerBBox(cvf::Vec3d(25.0, 25.0, -10.0), cvf::Vec3d(75.0, 75.0, 25.0));
+    overlapVertices = RigCellGeometryTools::estimateHexOverlapWithBoundingBox(cornerVertices, smallerInnerBBox, &overlapBoundingBox);
+    EXPECT_DOUBLE_EQ(50 * 50 * 25, RigCellGeometryTools::calculateCellVolume(overlapVertices));
+
+    cvf::BoundingBox smallerBBox(cvf::Vec3d(50.0, 50.0, 0.0), cvf::Vec3d(100.0, 100.0, 100.0));
+    overlapVertices = RigCellGeometryTools::estimateHexOverlapWithBoundingBox(cornerVertices, smallerBBox, &overlapBoundingBox);
+    double tipVolume = 50 * 50 * 50 * 0.5;
+    EXPECT_DOUBLE_EQ(tipVolume, RigCellGeometryTools::calculateCellVolume(overlapVertices));
+
+    cvf::BoundingBox smallerBBox2(cvf::Vec3d(0.0, 0.0, 0.0), cvf::Vec3d(50.0, 50.0, 100.0));
+    overlapVertices = RigCellGeometryTools::estimateHexOverlapWithBoundingBox(cornerVertices, smallerBBox2, &overlapBoundingBox);
+    double expectedVolume = (totalCellVolume - 2*tipVolume) * 0.5;
+    EXPECT_DOUBLE_EQ(expectedVolume, RigCellGeometryTools::calculateCellVolume(overlapVertices));
 }
 
 //--------------------------------------------------------------------------------------------------