///////////////////////////////////////////////////////////////////////////////// // // 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 // for more details. // ///////////////////////////////////////////////////////////////////////////////// #include "RigHexIntersectionTools.h" #include "../ReservoirDataModel/RigCellGeometryTools.h" #include "../LibGeometry/cvfBoundingBox.h" #include "../ReservoirDataModel/cvfGeometryTools.h" #include "../LibGeometry/cvfRay.h" #include "../cafHexGridIntersectionTools/cafHexGridIntersectionTools.h" namespace external { //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- int RigHexIntersectionTools::lineHexCellIntersection( const cvf::Vec3d p1, const cvf::Vec3d p2, const cvf::Vec3d hexCorners[8], const size_t hexIndex, std::vector* intersections ) { CVF_ASSERT( intersections != nullptr ); std::set uniqueIntersections; for ( int face = 0; face < 6; ++face ) { cvf::ubyte faceVertexIndices[4]; cvf::StructGridInterface::cellFaceVertexIndices( static_cast( face ), faceVertexIndices ); cvf::Vec3d intersection; bool isEntering = false; cvf::Vec3d faceCenter = cvf::GeometryTools::computeFaceCenter( hexCorners[faceVertexIndices[0]], hexCorners[faceVertexIndices[1]], hexCorners[faceVertexIndices[2]], hexCorners[faceVertexIndices[3]] ); for ( int i = 0; i < 4; ++i ) { int next = i < 3 ? i + 1 : 0; int intsStatus = cvf::GeometryTools::intersectLineSegmentTriangle( p1, p2, hexCorners[faceVertexIndices[i]], hexCorners[faceVertexIndices[next]], faceCenter, &intersection, &isEntering ); if ( intsStatus == 1 ) { uniqueIntersections.insert( HexIntersectionInfo( intersection, isEntering, static_cast( face ), hexIndex ) ); } } } int intersectionCount = 0; for ( const auto& intersection : uniqueIntersections ) { intersections->push_back( intersection ); ++intersectionCount; } return intersectionCount; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- bool RigHexIntersectionTools::lineIntersectsHexCell( const cvf::Vec3d p1, const cvf::Vec3d p2, const cvf::Vec3d hexCorners[8] ) { for ( int face = 0; face < 6; ++face ) { cvf::ubyte faceVertexIndices[4]; cvf::StructGridInterface::cellFaceVertexIndices( static_cast( face ), faceVertexIndices ); cvf::Vec3d intersection; bool isEntering = false; cvf::Vec3d faceCenter = cvf::GeometryTools::computeFaceCenter( hexCorners[faceVertexIndices[0]], hexCorners[faceVertexIndices[1]], hexCorners[faceVertexIndices[2]], hexCorners[faceVertexIndices[3]] ); for ( int i = 0; i < 4; ++i ) { int next = i < 3 ? i + 1 : 0; int intsStatus = cvf::GeometryTools::intersectLineSegmentTriangle( p1, p2, hexCorners[faceVertexIndices[i]], hexCorners[faceVertexIndices[next]], faceCenter, &intersection, &isEntering ); if ( intsStatus == 1 ) { return true; } } } return false; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- bool RigHexIntersectionTools::isPointInCell( const cvf::Vec3d point, const cvf::Vec3d hexCorners[8] ) { cvf::Ray ray; ray.setOrigin( point ); size_t intersections = 0; for ( int face = 0; face < 6; ++face ) { cvf::ubyte faceVertexIndices[4]; cvf::StructGridInterface::cellFaceVertexIndices( static_cast( face ), faceVertexIndices ); cvf::Vec3d faceCenter = cvf::GeometryTools::computeFaceCenter( hexCorners[faceVertexIndices[0]], hexCorners[faceVertexIndices[1]], hexCorners[faceVertexIndices[2]], hexCorners[faceVertexIndices[3]] ); for ( int i = 0; i < 4; ++i ) { int next = i < 3 ? i + 1 : 0; if ( ray.triangleIntersect( hexCorners[faceVertexIndices[i]], hexCorners[faceVertexIndices[next]], faceCenter ) ) { ++intersections; } } } return intersections % 2 == 1; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- // bool RigHexIntersectionTools::planeHexCellIntersection( cvf::Vec3d* hexCorners, // const cvf::Plane& fracturePlane, // std::list>& intersectionLineSegments ) // { // bool isCellIntersected = false; // cvf::ubyte faceVertexIndices[4]; // caf::HexGridIntersectionTools::ClipVx triangleIntersectionPoint1; // caf::HexGridIntersectionTools::ClipVx triangleIntersectionPoint2; // bool isMostVxesOnPositiveSideOfP1 = false; // for ( int face = 0; face < 6; ++face ) // { // cvf::StructGridInterface::cellFaceVertexIndices( static_cast( face ), // faceVertexIndices ); // cvf::Vec3d faceCenter = cvf::GeometryTools::computeFaceCenter( hexCorners[faceVertexIndices[0]], // hexCorners[faceVertexIndices[1]], // hexCorners[faceVertexIndices[2]], // hexCorners[faceVertexIndices[3]] ); // for ( int i = 0; i < 4; i++ ) // { // int next = i < 3 ? i + 1 : 0; // bool isIntersectingPlane = // caf::HexGridIntersectionTools::planeTriangleIntersection( fracturePlane, // hexCorners[faceVertexIndices[i]], // 0, // hexCorners[faceVertexIndices[next]], // 1, // faceCenter, // 2, // &triangleIntersectionPoint1, // &triangleIntersectionPoint2, // &isMostVxesOnPositiveSideOfP1 ); // if ( isIntersectingPlane ) // { // isCellIntersected = true; // intersectionLineSegments.emplace_back( triangleIntersectionPoint1.vx, triangleIntersectionPoint2.vx ); // } // } // } // return isCellIntersected; // } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- // bool RigHexIntersectionTools::planeHexIntersectionPolygons( std::array hexCorners, // cvf::Mat4d transformMatrixForPlane, // std::vector>& polygons ) // { // bool isCellIntersected = false; // cvf::Plane fracturePlane; // fracturePlane.setFromPointAndNormal( transformMatrixForPlane.translation(), // static_cast( transformMatrixForPlane.col( 2 ) ) ); // // Find line-segments where cell and fracture plane intersects // std::list> intersectionLineSegments; // isCellIntersected = planeHexCellIntersection( hexCorners.data(), fracturePlane, intersectionLineSegments ); // RigCellGeometryTools::createPolygonFromLineSegments( intersectionLineSegments, polygons ); // return isCellIntersected; // } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- bool operator<( const HexIntersectionInfo& hi1, const HexIntersectionInfo& hi2 ) { const double tolerance = 1e-6; if ( hi1.m_hexIndex != hi2.m_hexIndex ) return hi1.m_hexIndex < hi2.m_hexIndex; if ( hi1.m_face != hi2.m_face ) return hi1.m_face < hi2.m_face; if ( hi1.m_isIntersectionEntering != hi2.m_isIntersectionEntering ) return hi1.m_isIntersectionEntering < hi2.m_isIntersectionEntering; if ( hi1.m_face != hi2.m_face ) return hi1.m_face < hi2.m_face; if ( cvf::Math::abs( hi2.m_intersectionPoint.x() - hi1.m_intersectionPoint.x() ) > tolerance ) return hi1.m_intersectionPoint.x() < hi2.m_intersectionPoint.x(); if ( cvf::Math::abs( hi2.m_intersectionPoint.y() - hi1.m_intersectionPoint.y() ) > tolerance ) return hi1.m_intersectionPoint.y() < hi2.m_intersectionPoint.y(); if ( cvf::Math::abs( hi2.m_intersectionPoint.z() - hi1.m_intersectionPoint.z() ) > tolerance ) return hi1.m_intersectionPoint.z() < hi2.m_intersectionPoint.z(); return false; } } //namespace external