#1091 - pre-proto - Moving functions planeHexCellIntersection and createPolygonFromLineSegments from RimFracture to RigCellGeometryTools

This commit is contained in:
astridkbjorke 2017-01-23 15:22:51 +01:00
parent eb423dd879
commit 02f4b94ca1
5 changed files with 188 additions and 104 deletions

View File

@ -21,6 +21,7 @@
#include "RiaApplication.h"
#include "RigCell.h"
#include "RigCellGeometryTools.h"
#include "RigFracture.h"
#include "RigMainGrid.h"
#include "RigTesselatorTools.h"
@ -273,7 +274,7 @@ void RimFracture::computeTransmissibility()
//Transform planCell polygon(s) to x/y coordinate system (where fracturePolygon already is located)
cvf::Mat4f invertedTransMatrix = transformMatrix().getInverted();
for (std::vector<cvf::Vec3d> planeCellPolygon : planeCellPolygons)
for (std::vector<cvf::Vec3d> & planeCellPolygon : planeCellPolygons)
{
for (cvf::Vec3d& v : planeCellPolygon)
{
@ -367,110 +368,13 @@ bool RimFracture::planeCellIntersectionPolygons(size_t cellindex, std::vector<st
//Find line-segments where cell and fracture plane intersects
std::list<std::pair<cvf::Vec3d, cvf::Vec3d > > intersectionLineSegments;
isCellIntersected = planeHexCellIntersection(hexCorners, fracturePlane, intersectionLineSegments);
isCellIntersected = RigCellGeometryTools::planeHexCellIntersection(hexCorners, fracturePlane, intersectionLineSegments);
createPolygonFromLineSegments(intersectionLineSegments, polygons);
RigCellGeometryTools::createPolygonFromLineSegments(intersectionLineSegments, polygons);
return isCellIntersected;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RimFracture::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;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RimFracture::createPolygonFromLineSegments(std::list<std::pair<cvf::Vec3d, cvf::Vec3d>> &intersectionLineSegments, std::vector<std::vector<cvf::Vec3d>> &polygons)
{
bool startNewPolygon = true;
while (!intersectionLineSegments.empty())
{
if (startNewPolygon)
{
std::vector<cvf::Vec3d> polygon;
//Add first line segments to polygon and remove from list
std::pair<cvf::Vec3d, cvf::Vec3d > linesegment = intersectionLineSegments.front();
polygon.push_back(linesegment.first);
polygon.push_back(linesegment.second);
intersectionLineSegments.remove(linesegment);
polygons.push_back(polygon);
startNewPolygon = false;
}
std::vector<cvf::Vec3d>& polygon = polygons.back();
//Search remaining list for next point...
bool isFound = false;
float tolerance = 0.0001f;
for (std::list<std::pair<cvf::Vec3d, cvf::Vec3d > >::iterator lIt = intersectionLineSegments.begin(); lIt != intersectionLineSegments.end(); lIt++)
{
cvf::Vec3d lineSegmentStart = lIt->first;
cvf::Vec3d lineSegmentEnd = lIt->second;
cvf::Vec3d polygonEnd = polygon.back();
if (((lineSegmentStart - polygonEnd).lengthSquared() < tolerance*tolerance))
{
polygon.push_back(lIt->second);
intersectionLineSegments.erase(lIt);
isFound = true;
break;
}
if (((lineSegmentEnd - polygonEnd).lengthSquared() < tolerance*tolerance))
{
polygon.push_back(lIt->first);
intersectionLineSegments.erase(lIt);
isFound = true;
break;
}
}
if (isFound)
{
continue;
}
else
{
startNewPolygon = true;
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

View File

@ -94,10 +94,6 @@ private:
//Functions for area calculations - should these be in separate class
bool planeCellIntersectionPolygons(size_t cellindex, std::vector<std::vector<cvf::Vec3d> > & polygons);
bool planeHexCellIntersection(cvf::Vec3d * hexCorners, cvf::Plane fracturePlane, std::list<std::pair<cvf::Vec3d, cvf::Vec3d > > & intersectionLineSegments);
void createPolygonFromLineSegments(std::list<std::pair<cvf::Vec3d, cvf::Vec3d>> &intersectionLineSegments, std::vector<std::vector<cvf::Vec3d>> &polygons);
void calculateFracturePlaneCellPolygonOverlapArea(std::vector<std::vector<cvf::Vec3d> > planeCellPolygons,
std::vector<cvf::Vec3f> fracturePolygon, double & area);

View File

@ -51,6 +51,7 @@ ${CEE_CURRENT_LIST_DIR}RigSummaryCaseData.h
${CEE_CURRENT_LIST_DIR}RigLasFileExporter.h
${CEE_CURRENT_LIST_DIR}RigFracture.h
${CEE_CURRENT_LIST_DIR}RigTesselatorTools.h
${CEE_CURRENT_LIST_DIR}RigCellGeometryTools.h
)
set (SOURCE_GROUP_SOURCE_FILES
@ -96,6 +97,8 @@ ${CEE_CURRENT_LIST_DIR}RigSummaryCaseData.cpp
${CEE_CURRENT_LIST_DIR}RigLasFileExporter.cpp
${CEE_CURRENT_LIST_DIR}RigFracture.cpp
${CEE_CURRENT_LIST_DIR}RigTesselatorTools.cpp
${CEE_CURRENT_LIST_DIR}RigCellGeometryTools.cpp
)
list(APPEND CODE_HEADER_FILES

View File

@ -0,0 +1,141 @@
/////////////////////////////////////////////////////////////////////////////////
//
// 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 "RigCellGeometryTools.h"
#include "cvfStructGrid.h"
#include "cvfGeometryTools.h"
#include "cafHexGridIntersectionTools\cafHexGridIntersectionTools.h"
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigCellGeometryTools::RigCellGeometryTools()
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigCellGeometryTools::~RigCellGeometryTools()
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
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;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigCellGeometryTools::createPolygonFromLineSegments(std::list<std::pair<cvf::Vec3d, cvf::Vec3d>> &intersectionLineSegments, std::vector<std::vector<cvf::Vec3d>> &polygons)
{
bool startNewPolygon = true;
while (!intersectionLineSegments.empty())
{
if (startNewPolygon)
{
std::vector<cvf::Vec3d> polygon;
//Add first line segments to polygon and remove from list
std::pair<cvf::Vec3d, cvf::Vec3d > linesegment = intersectionLineSegments.front();
polygon.push_back(linesegment.first);
polygon.push_back(linesegment.second);
intersectionLineSegments.remove(linesegment);
polygons.push_back(polygon);
startNewPolygon = false;
}
std::vector<cvf::Vec3d>& polygon = polygons.back();
//Search remaining list for next point...
bool isFound = false;
float tolerance = 0.0001f;
for (std::list<std::pair<cvf::Vec3d, cvf::Vec3d > >::iterator lIt = intersectionLineSegments.begin(); lIt != intersectionLineSegments.end(); lIt++)
{
cvf::Vec3d lineSegmentStart = lIt->first;
cvf::Vec3d lineSegmentEnd = lIt->second;
cvf::Vec3d polygonEnd = polygon.back();
if (((lineSegmentStart - polygonEnd).lengthSquared() < tolerance*tolerance))
{
polygon.push_back(lIt->second);
intersectionLineSegments.erase(lIt);
isFound = true;
break;
}
if (((lineSegmentEnd - polygonEnd).lengthSquared() < tolerance*tolerance))
{
polygon.push_back(lIt->first);
intersectionLineSegments.erase(lIt);
isFound = true;
break;
}
}
if (isFound)
{
continue;
}
else
{
startNewPolygon = true;
}
}
}

View File

@ -0,0 +1,40 @@
/////////////////////////////////////////////////////////////////////////////////
//
// 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.
//
/////////////////////////////////////////////////////////////////////////////////
#pragma once
#include "cvfBase.h"
#include "cvfVector3.h"
#include "cvfPlane.h"
#include <list>
class RigCellGeometryTools
{
public:
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);
};