diff --git a/ApplicationCode/ProjectDataModel/Rim2dGridProjection.cpp b/ApplicationCode/ProjectDataModel/Rim2dGridProjection.cpp index f954560fa2..b85424b788 100644 --- a/ApplicationCode/ProjectDataModel/Rim2dGridProjection.cpp +++ b/ApplicationCode/ProjectDataModel/Rim2dGridProjection.cpp @@ -684,7 +684,7 @@ std::vector> Rim2dGridProjection::visibleCellsAndWeigh localGrid->cellCornerVertices(localCellIdx, hexCorners.data()); cvf::BoundingBox overlapBBox; - std::array overlapCorners = createHexOverlapEstimation(bbox2dElement, hexCorners, &overlapBBox); + std::array 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 Rim2dGridProjection::createHexOverlapEstimation(const cvf::BoundingBox& bbox2dElement, const std::array& hexCorners, cvf::BoundingBox* overlapBoundingBox) const -{ - CVF_ASSERT(overlapBoundingBox); - *overlapBoundingBox = cvf::BoundingBox(); - std::array 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; -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp b/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp index 34c52fafcc..6bf94fe9fa 100644 --- a/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp +++ b/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp @@ -27,6 +27,8 @@ #include "clipper/clipper.hpp" +#include "cvfMath.h" + #include #include @@ -82,6 +84,43 @@ double RigCellGeometryTools::calculateCellVolume(const std::array return std::abs(volume); // Altogether 18 + 3*17 + 3 + 1 flops = 73 flops. } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::array RigCellGeometryTools::estimateHexOverlapWithBoundingBox(const std::array& hexCorners, const cvf::BoundingBox& boundingBox, cvf::BoundingBox* overlapBoundingBox) +{ + CVF_ASSERT(overlapBoundingBox); + *overlapBoundingBox = cvf::BoundingBox(); + std::array 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; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h b/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h index f54b4d6969..9de7ae10d8 100644 --- a/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h +++ b/ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h @@ -1,17 +1,17 @@ ///////////////////////////////////////////////////////////////////////////////// // // Copyright (C) 2017 Statoil ASA -// +// // ResInsight is free software: you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation, either version 3 of the License, or // (at your option) any later version. -// +// // ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY // WARRANTY; without even the implied warranty of MERCHANTABILITY or // FITNESS FOR A PARTICULAR PURPOSE. -// -// See the GNU General Public License at +// +// See the GNU General Public License at // for more details. // ///////////////////////////////////////////////////////////////////////////////// @@ -28,38 +28,50 @@ #include #include - - class RigCellGeometryTools { public: - static double calculateCellVolume(const std::array& hexCorners); + static double calculateCellVolume(const std::array& hexCorners); + static std::array estimateHexOverlapWithBoundingBox(const std::array& hexCorners, + const cvf::BoundingBox& boundingBox2dExtrusion, + cvf::BoundingBox* overlapBoundingBox); - static void createPolygonFromLineSegments(std::list> &intersectionLineSegments, std::vector> &polygons); + static void createPolygonFromLineSegments(std::list>& intersectionLineSegments, + std::vector>& polygons); - static void findCellLocalXYZ(const std::array& hexCorners, cvf::Vec3d &localXdirection, cvf::Vec3d &localYdirection, cvf::Vec3d &localZdirection); + static void findCellLocalXYZ(const std::array& hexCorners, + cvf::Vec3d& localXdirection, + cvf::Vec3d& localYdirection, + cvf::Vec3d& localZdirection); static double polygonLengthInLocalXdirWeightedByArea(const std::vector& polygon2d); static std::vector> intersectPolygons(const std::vector& polygon1, const std::vector& polygon2); - + static std::vector> subtractPolygons(const std::vector& sourcePolygon, const std::vector>& polygonToSubtract); - enum ZInterpolationType { INTERPOLATE_LINE_Z, USE_HUGEVAL, USE_ZERO}; - static std::vector > clipPolylineByPolygon(const std::vector& polyLine, - const std::vector& polygon, - ZInterpolationType interpolType = USE_ZERO); + enum ZInterpolationType + { + INTERPOLATE_LINE_Z, + USE_HUGEVAL, + USE_ZERO + }; + static std::vector> clipPolylineByPolygon(const std::vector& polyLine, + const std::vector& polygon, + ZInterpolationType interpolType = USE_ZERO); - static std::pair getLineThroughBoundingBox(const cvf::Vec3d& lineDirection, const cvf::BoundingBox& polygonBBox, const cvf::Vec3d& pointOnLine); + static std::pair getLineThroughBoundingBox(const cvf::Vec3d& lineDirection, + const cvf::BoundingBox& polygonBBox, + const cvf::Vec3d& pointOnLine); - static double getLengthOfPolygonAlongLine(const std::pair& line, const std::vector& polygon); + static double getLengthOfPolygonAlongLine(const std::pair& line, + const std::vector& polygon); static std::vector unionOfPolygons(const std::vector>& polygons); private: - static std::vector ajustPolygonToAvoidIntersectionsAtVertex(const std::vector& polyLine, + static std::vector ajustPolygonToAvoidIntersectionsAtVertex(const std::vector& polyLine, const std::vector& polygon); - }; diff --git a/ApplicationCode/UnitTests/RigCellGeometryTools-Test.cpp b/ApplicationCode/UnitTests/RigCellGeometryTools-Test.cpp index 366abf65c5..3ec9fbdc24 100644 --- a/ApplicationCode/UnitTests/RigCellGeometryTools-Test.cpp +++ b/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 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; - - double extraVolume = 0.5 * extent.z() * 100 * extent.y(); + // 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() * 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 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 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 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 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)); } //--------------------------------------------------------------------------------------------------