Skip to content

Commit

Permalink
#1602 Consolidate some intersection functionality to common tools class
Browse files Browse the repository at this point in the history
  • Loading branch information
bjornarg committed Jun 16, 2017
1 parent 2e764b3 commit 2090db6
Show file tree
Hide file tree
Showing 9 changed files with 125 additions and 99 deletions.
36 changes: 0 additions & 36 deletions ApplicationCode/ReservoirDataModel/RigCellGeometryTools.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -29,42 +29,6 @@
#include <vector>


//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigCellGeometryTools::planeHexCellIntersection(cvf::Vec3d * hexCorners, cvf::Plane fracturePlane, std::list<std::pair<cvf::Vec3d, cvf::Vec3d > > & intersectionLineSegments)
{
bool isCellIntersected = false;
for (int face = 0; face < 6; ++face)
{
cvf::ubyte faceVertexIndices[4];
cvf::StructGridInterface::cellFaceVertexIndices(static_cast<cvf::StructGridInterface::FaceType>(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;
caf::HexGridIntersectionTools::ClipVx triangleIntersectionPoint1;
caf::HexGridIntersectionTools::ClipVx triangleIntersectionPoint2;

bool isMostVxesOnPositiveSideOfP1 = false;

bool isIntersectingPlane = caf::HexGridIntersectionTools::planeTriangleIntersection(fracturePlane,
hexCorners[faceVertexIndices[i]], 0,
hexCorners[faceVertexIndices[next]], 1,
faceCenter, 2,
&triangleIntersectionPoint1, &triangleIntersectionPoint2, &isMostVxesOnPositiveSideOfP1);

if (isIntersectingPlane)
{
isCellIntersected = true;
intersectionLineSegments.push_back({ triangleIntersectionPoint1.vx, triangleIntersectionPoint2.vx });
}
}
} return isCellIntersected;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
Expand Down
2 changes: 0 additions & 2 deletions ApplicationCode/ReservoirDataModel/RigCellGeometryTools.h
Original file line number Diff line number Diff line change
Expand Up @@ -36,8 +36,6 @@ class RigCellGeometryTools
RigCellGeometryTools();
~RigCellGeometryTools();

static bool planeHexCellIntersection(cvf::Vec3d * hexCorners, cvf::Plane fracturePlane, std::list<std::pair<cvf::Vec3d, cvf::Vec3d > > & intersectionLineSegments);

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);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
#include "RigFractureTransmissibilityEquations.h"
#include "RigMainGrid.h"
#include "RigResultAccessorFactory.h"
#include "RigHexIntersectionTools.h"

#include "RimEclipseCase.h"
#include "RimReservoirCellResultsStorage.h"
Expand Down Expand Up @@ -134,7 +135,7 @@ void RigEclipseToStimPlanCellTransmissibilityCalculator::calculateStimPlanCellsM
mainGrid->cellCornerVertices(fracCell, hexCorners.data());

std::vector<std::vector<cvf::Vec3d> > planeCellPolygons;
bool isPlanIntersected = planeCellIntersectionPolygons(hexCorners.data(), m_fractureTransform, planeCellPolygons);
bool isPlanIntersected = RigHexIntersectionTools::planeHexIntersectionPolygons(hexCorners, m_fractureTransform, planeCellPolygons);
if (!isPlanIntersected || planeCellPolygons.size() == 0) continue;

cvf::Vec3d localX;
Expand Down Expand Up @@ -241,28 +242,3 @@ std::vector<size_t> RigEclipseToStimPlanCellTransmissibilityCalculator::getPoten

return cellIndices;
}


//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigEclipseToStimPlanCellTransmissibilityCalculator::planeCellIntersectionPolygons(cvf::Vec3d hexCorners[8],
cvf::Mat4f transformMatrixForPlane,
std::vector<std::vector<cvf::Vec3d> > & polygons)
{
bool isCellIntersected = false;

cvf::Plane fracturePlane;
fracturePlane.setFromPointAndNormal(static_cast<cvf::Vec3d>(transformMatrixForPlane.translation()),
static_cast<cvf::Vec3d>(transformMatrixForPlane.col(2)));

//Find line-segments where cell and fracture plane intersects
std::list<std::pair<cvf::Vec3d, cvf::Vec3d > > intersectionLineSegments;

isCellIntersected = RigCellGeometryTools::planeHexCellIntersection(hexCorners, fracturePlane, intersectionLineSegments);

RigCellGeometryTools::createPolygonFromLineSegments(intersectionLineSegments, polygons);

return isCellIntersected;
}

Original file line number Diff line number Diff line change
Expand Up @@ -47,10 +47,6 @@ class RigEclipseToStimPlanCellTransmissibilityCalculator
std::vector<size_t> getPotentiallyFracturedCellsForPolygon(std::vector<cvf::Vec3d> polygon);


static bool planeCellIntersectionPolygons(cvf::Vec3d hexCorners[8],
cvf::Mat4f transformMatrixForPlane,
std::vector<std::vector<cvf::Vec3d> > & polygons);

private:
const RimEclipseCase* m_case;
double m_cDarcy;
Expand Down
63 changes: 63 additions & 0 deletions ApplicationCode/ReservoirDataModel/RigHexIntersectionTools.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -17,10 +17,14 @@
/////////////////////////////////////////////////////////////////////////////////
#include "RigHexIntersectionTools.h"

#include "RigCellGeometryTools.h"

#include "cvfBoundingBox.h"
#include "cvfGeometryTools.h"
#include "cvfRay.h"

#include "cafHexGridIntersectionTools/cafHexGridIntersectionTools.h"


//--------------------------------------------------------------------------------------------------
///
Expand Down Expand Up @@ -89,3 +93,62 @@ bool RigHexIntersectionTools::isPointInCell(const cvf::Vec3d point, const cvf::V
}
return intersections % 2 == 1;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigHexIntersectionTools::planeHexCellIntersection(cvf::Vec3d* hexCorners, cvf::Plane fracturePlane, std::list<std::pair<cvf::Vec3d, cvf::Vec3d > >& intersectionLineSegments)
{
bool isCellIntersected = false;
for (int face = 0; face < 6; ++face)
{
cvf::ubyte faceVertexIndices[4];
cvf::StructGridInterface::cellFaceVertexIndices(static_cast<cvf::StructGridInterface::FaceType>(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;
caf::HexGridIntersectionTools::ClipVx triangleIntersectionPoint1;
caf::HexGridIntersectionTools::ClipVx triangleIntersectionPoint2;

bool isMostVxesOnPositiveSideOfP1 = false;

bool isIntersectingPlane = caf::HexGridIntersectionTools::planeTriangleIntersection(fracturePlane,
hexCorners[faceVertexIndices[i]], 0,
hexCorners[faceVertexIndices[next]], 1,
faceCenter, 2,
&triangleIntersectionPoint1, &triangleIntersectionPoint2, &isMostVxesOnPositiveSideOfP1);

if (isIntersectingPlane)
{
isCellIntersected = true;
intersectionLineSegments.push_back({ triangleIntersectionPoint1.vx, triangleIntersectionPoint2.vx });
}
}
} return isCellIntersected;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigHexIntersectionTools::planeHexIntersectionPolygons(std::array<cvf::Vec3d, 8> hexCorners,
cvf::Mat4f transformMatrixForPlane,
std::vector<std::vector<cvf::Vec3d> >& polygons)
{
bool isCellIntersected = false;

cvf::Plane fracturePlane;
fracturePlane.setFromPointAndNormal(static_cast<cvf::Vec3d>(transformMatrixForPlane.translation()),
static_cast<cvf::Vec3d>(transformMatrixForPlane.col(2)));

//Find line-segments where cell and fracture plane intersects
std::list<std::pair<cvf::Vec3d, cvf::Vec3d > > intersectionLineSegments;

isCellIntersected = planeHexCellIntersection(hexCorners.data(), fracturePlane, intersectionLineSegments);

RigCellGeometryTools::createPolygonFromLineSegments(intersectionLineSegments, polygons);

return isCellIntersected;
}
9 changes: 9 additions & 0 deletions ApplicationCode/ReservoirDataModel/RigHexIntersectionTools.h
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,9 @@
#include "cvfBase.h"
#include "cvfVector3.h"
#include "cvfStructGrid.h"
#include "cvfPlane.h"

#include <array>

//==================================================================================================
/// Internal class for intersection point info
Expand Down Expand Up @@ -58,7 +61,13 @@ struct RigHexIntersectionTools

static bool isPointInCell(const cvf::Vec3d point, const cvf::Vec3d hexCorners[8]);

static bool planeHexCellIntersection(cvf::Vec3d* hexCorners,
cvf::Plane fracturePlane,
std::list<std::pair<cvf::Vec3d, cvf::Vec3d > >& intersectionLineSegments);

static bool planeHexIntersectionPolygons(std::array<cvf::Vec3d, 8> hexCorners,
cvf::Mat4f transformMatrixForPlane,
std::vector<std::vector<cvf::Vec3d> > & polygons);
};


1 change: 1 addition & 0 deletions ApplicationCode/UnitTests/CMakeLists_files.cmake
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@ ${CEE_CURRENT_LIST_DIR}WellPathAsciiFileReader-Test.cpp
${CEE_CURRENT_LIST_DIR}opm-flowdiagnostics-Test.cpp
${CEE_CURRENT_LIST_DIR}RigCellGeometryTools-Test.cpp
${CEE_CURRENT_LIST_DIR}RigTransmissibilityCondenser-Test.cpp
${CEE_CURRENT_LIST_DIR}RigHexIntersectionTools-Test.cpp

)

Expand Down
31 changes: 0 additions & 31 deletions ApplicationCode/UnitTests/RigCellGeometryTools-Test.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -104,37 +104,6 @@ TEST(RigCellGeometryTools, createPolygonTestRealCase)

}



//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
TEST(RigCellGeometryTools, planeHexCellIntersectionTest)
{
cvf::Vec3d hexCorners[8];
hexCorners[0] = cvf::Vec3d(0, 0, 0);
hexCorners[1] = cvf::Vec3d(1, 0, 0);
hexCorners[2] = cvf::Vec3d(0, 1, 0);
hexCorners[3] = cvf::Vec3d(0, 0, 1);
hexCorners[4] = cvf::Vec3d(0, 1, 1);
hexCorners[5] = cvf::Vec3d(1, 1, 0);
hexCorners[6] = cvf::Vec3d(1, 0, 1);
hexCorners[7] = cvf::Vec3d(1, 1, 1);

std::list<std::pair<cvf::Vec3d, cvf::Vec3d > > intersectionLineSegments;
bool isCellIntersected = false;
cvf::Plane fracturePlane;

fracturePlane.setFromPointAndNormal(cvf::Vec3d(0.5, 0.5, 0.5), cvf::Vec3d(1, 0, 0) );
isCellIntersected = RigCellGeometryTools::planeHexCellIntersection(hexCorners, fracturePlane, intersectionLineSegments);
EXPECT_TRUE(isCellIntersected);

fracturePlane.setFromPointAndNormal(cvf::Vec3d(1.5, 1.5, 1.5), cvf::Vec3d(1, 0, 0));
isCellIntersected = RigCellGeometryTools::planeHexCellIntersection(hexCorners, fracturePlane, intersectionLineSegments);
EXPECT_FALSE(isCellIntersected);

}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
Expand Down
50 changes: 50 additions & 0 deletions ApplicationCode/UnitTests/RigHexIntersectionTools-Test.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@
/////////////////////////////////////////////////////////////////////////////////
//
// 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 <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////

#include "gtest/gtest.h"

#include "RigHexIntersectionTools.h"

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
TEST(RigHexIntersectionTools, planeHexCellIntersectionTest)
{
cvf::Vec3d hexCorners[8];
hexCorners[0] = cvf::Vec3d(0, 0, 0);
hexCorners[1] = cvf::Vec3d(1, 0, 0);
hexCorners[2] = cvf::Vec3d(0, 1, 0);
hexCorners[3] = cvf::Vec3d(0, 0, 1);
hexCorners[4] = cvf::Vec3d(0, 1, 1);
hexCorners[5] = cvf::Vec3d(1, 1, 0);
hexCorners[6] = cvf::Vec3d(1, 0, 1);
hexCorners[7] = cvf::Vec3d(1, 1, 1);

std::list<std::pair<cvf::Vec3d, cvf::Vec3d > > intersectionLineSegments;
bool isCellIntersected = false;
cvf::Plane fracturePlane;

fracturePlane.setFromPointAndNormal(cvf::Vec3d(0.5, 0.5, 0.5), cvf::Vec3d(1, 0, 0) );
isCellIntersected = RigHexIntersectionTools::planeHexCellIntersection(hexCorners, fracturePlane, intersectionLineSegments);
EXPECT_TRUE(isCellIntersected);

fracturePlane.setFromPointAndNormal(cvf::Vec3d(1.5, 1.5, 1.5), cvf::Vec3d(1, 0, 0));
isCellIntersected = RigHexIntersectionTools::planeHexCellIntersection(hexCorners, fracturePlane, intersectionLineSegments);
EXPECT_FALSE(isCellIntersected);

}

0 comments on commit 2090db6

Please sign in to comment.