From e171aa66acd0fec857c8798abf5485050cdf1c71 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jacob=20St=C3=B8ren?= Date: Wed, 11 Dec 2013 14:55:14 +0100 Subject: [PATCH] NNC processing first cut --- .../FileInterface/RifReaderEclipseOutput.cpp | 19 ++- .../ReservoirDataModel/RigNNCData.cpp | 151 +++++++++++------- .../ReservoirDataModel/RigNNCData.h | 15 +- 3 files changed, 110 insertions(+), 75 deletions(-) diff --git a/ApplicationCode/FileInterface/RifReaderEclipseOutput.cpp b/ApplicationCode/FileInterface/RifReaderEclipseOutput.cpp index 0d90d8c556..5c3c4c1f5d 100644 --- a/ApplicationCode/FileInterface/RifReaderEclipseOutput.cpp +++ b/ApplicationCode/FileInterface/RifReaderEclipseOutput.cpp @@ -392,15 +392,22 @@ bool RifReaderEclipseOutput::open(const QString& fileName, RigCaseData* eclipseC m_eclipseCase = eclipseCase; - progInfo.setProgressDescription("Reading Result index"); - progInfo.setNextProgressIncrement(50); - // Build results meta data + progInfo.setProgressDescription("Reading Result index"); + progInfo.setNextProgressIncrement(25); buildMetaData(); - - transferNNCData(mainEclGrid, m_ecl_init_file, eclipseCase->mainGrid()); - progInfo.incrementProgress(); + + progInfo.setProgressDescription("Reading NNC data"); + progInfo.setNextProgressIncrement(5); + transferNNCData(mainEclGrid, m_ecl_init_file, eclipseCase->mainGrid()); + progInfo.incrementProgress(); + + progInfo.setProgressDescription("Processing NNC's"); + progInfo.setNextProgressIncrement(20); + eclipseCase->mainGrid()->nncData()->processConnections( *(eclipseCase->mainGrid())); + progInfo.incrementProgress(); + progInfo.setNextProgressIncrement(8); progInfo.setProgressDescription("Reading Well information"); readWellCells(mainEclGrid); diff --git a/ApplicationCode/ReservoirDataModel/RigNNCData.cpp b/ApplicationCode/ReservoirDataModel/RigNNCData.cpp index 8d81047055..8081ba644b 100644 --- a/ApplicationCode/ReservoirDataModel/RigNNCData.cpp +++ b/ApplicationCode/ReservoirDataModel/RigNNCData.cpp @@ -18,7 +18,7 @@ #include "RigNNCData.h" #include "RigMainGrid.h" -//#include "../ModelVisualization/cvfGeometryTools.h" +#include "cvfGeometryTools.h" //-------------------------------------------------------------------------------------------------- @@ -35,14 +35,16 @@ RigNNCData::RigNNCData() //-------------------------------------------------------------------------------------------------- void RigNNCData::processConnections(const RigMainGrid& mainGrid) { - /* - for (size_t cnIdx = 0; cnIdx < 0; ++cnIdx) + + for (size_t cnIdx = 0; cnIdx < m_connections.size(); ++cnIdx) { const RigCell& c1 = mainGrid.cells()[m_connections[cnIdx].m_c1GlobIdx]; const RigCell& c2 = mainGrid.cells()[m_connections[cnIdx].m_c2GlobIdx]; // Try to find the shared face + char hasNeighbourInAnyDirection = 0; + bool isPossibleNeighborInDirection[6]= {true, true, true, true, true, true}; if (c1.hostGrid() == c2.hostGrid()) { @@ -51,7 +53,7 @@ void RigNNCData::processConnections(const RigMainGrid& mainGrid) size_t i2, j2, k2; c2.hostGrid()->ijkFromCellIndex(c2.cellIndex(), &i2, &j2, &k2); - bool isPossibleNeighborInDirection[6]= {false, false, false, false, false, false}; + isPossibleNeighborInDirection[cvf::StructGridInterface::POS_I] = ((i1 + 1) == i2); isPossibleNeighborInDirection[cvf::StructGridInterface::NEG_I] = ((i2 + 1) == i1); isPossibleNeighborInDirection[cvf::StructGridInterface::POS_J] = ((j1 + 1) == j2); @@ -60,7 +62,7 @@ void RigNNCData::processConnections(const RigMainGrid& mainGrid) isPossibleNeighborInDirection[cvf::StructGridInterface::NEG_K] = ((k2 + 1) == k1); hasNeighbourInAnyDirection = - isPossibleNeighborInDirection[cvf::StructGridInterface::POS_I] + isPossibleNeighborInDirection[cvf::StructGridInterface::POS_I] + isPossibleNeighborInDirection[cvf::StructGridInterface::NEG_I] + isPossibleNeighborInDirection[cvf::StructGridInterface::POS_J] + isPossibleNeighborInDirection[cvf::StructGridInterface::NEG_J] @@ -73,68 +75,95 @@ void RigNNCData::processConnections(const RigMainGrid& mainGrid) if (!hasNeighbourInAnyDirection) { - m_connections[cnIdx].m_hasNoSharedArea = true; + // Add to search map + m_cellIdxToFaceToConnectionIdxMap[m_connections[cnIdx].m_c1GlobIdx][cvf::StructGridInterface::NO_FACE].push_back(cnIdx); + m_cellIdxToFaceToConnectionIdxMap[m_connections[cnIdx].m_c2GlobIdx][cvf::StructGridInterface::NO_FACE].push_back(cnIdx); + continue; // to next connection } + } - if (hasNeighbourInAnyDirection == 1) + // Possibly do some testing to avoid unneccesary overlap calculations + + cvf::Vec3d normal; + for (char fIdx = 0; fIdx < 6; ++fIdx) + { + if (isPossibleNeighborInDirection[fIdx]) { - for (char fIdx = 0; fIdx < 6; ++fIdx) + cvf::Vec3d fc1 = c1.faceCenter((cvf::StructGridInterface::FaceType)(fIdx)); + cvf::Vec3d fc2 = c2.faceCenter(cvf::StructGridInterface::oppositeFace((cvf::StructGridInterface::FaceType)(fIdx))); + cvf::Vec3d fc1ToFc2 = fc2 - fc1; + normal = c1.faceNormal((cvf::StructGridInterface::FaceType)(fIdx)); + normal.normalize(); + // Check that face centers are approx in the face plane + if (normal.dot(fc1ToFc2) < 0.01*fc1ToFc2.length()) { - if (isPossibleNeighborInDirection[fIdx]) - { - m_connections[cnIdx].m_c1Face = (cvf::StructGridInterface::FaceType)fIdx; - break; // the face loop - } + } - - // calculate polygon for this face - - } - else - { - - cvf::Vec3d normal; - for (char fIdx = 0; fIdx < 6; ++fIdx) - { - if (isPossibleNeighborInDirection[fIdx]) - { - cvf::Vec3d fc1 = c1.faceCenter((cvf::StructGridInterface::FaceType)(fIdx)); - cvf::Vec3d fc2 = c2.faceCenter(cvf::StructGridInterface::oppositeFace((cvf::StructGridInterface::FaceType)(fIdx))); - cvf::Vec3d fc1ToFc2 = fc2 - fc1; - normal = c1.faceNormal((cvf::StructGridInterface::FaceType)(fIdx)); - normal.normalize(); - // Check that face centers are approx in the face plane - if (normal.dot(fc1ToFc2) < 0.01*fc1ToFc2.length()) - { - - } - - // Calculate connection polygon - - std::vector polygon; - std::vector intersections; - bool isOk = false; - caf::SizeTArray4 face1; - caf::SizeTArray4 face2; - c1.faceIndices((cvf::StructGridInterface::FaceType)(fIdx), &face1); - c2.faceIndices(cvf::StructGridInterface::oppositeFace((cvf::StructGridInterface::FaceType)(fIdx)), &face2); - - isOk = cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads( - &polygon, - &intersections, - (cvf::EdgeIntersectStorage*)NULL, - cvf::wrapArrayConst(&mainGrid.nodes()), - face1.data(), - face2.data(), - 1e-6); - - } - - } - } - } } - */ + + for (char fIdx = 0; fIdx < 6; ++fIdx) + { + if (!isPossibleNeighborInDirection[fIdx]) + { + continue; + } + + // Calculate connection polygon + + std::vector polygon; + std::vector intersections; + caf::SizeTArray4 face1; + caf::SizeTArray4 face2; + c1.faceIndices((cvf::StructGridInterface::FaceType)(fIdx), &face1); + c2.faceIndices(cvf::StructGridInterface::oppositeFace((cvf::StructGridInterface::FaceType)(fIdx)), &face2); + + bool foundOverlap = cvf::GeometryTools::calculateOverlapPolygonOfTwoQuads( + &polygon, + &intersections, + (cvf::EdgeIntersectStorage*)NULL, + cvf::wrapArrayConst(&mainGrid.nodes()), + face1.data(), + face2.data(), + 1e-6); + + if (foundOverlap) + { + // Found an overlap polygon. Store data about connection + + m_connections[cnIdx].m_c1Face = (cvf::StructGridInterface::FaceType)fIdx; + for (size_t pIdx = 0; pIdx < polygon.size(); ++pIdx) + { + if (polygon[pIdx] < mainGrid.nodes().size()) + m_connections[cnIdx].m_polygon.push_back(mainGrid.nodes()[polygon[pIdx]]); + else + m_connections[cnIdx].m_polygon.push_back(intersections[polygon[pIdx] - mainGrid.nodes().size()]); + } + + // Add to search map + m_cellIdxToFaceToConnectionIdxMap[m_connections[cnIdx].m_c1GlobIdx][fIdx].push_back(cnIdx); + m_cellIdxToFaceToConnectionIdxMap[m_connections[cnIdx].m_c2GlobIdx][cvf::StructGridInterface::oppositeFace((cvf::StructGridInterface::FaceType)(fIdx))].push_back(cnIdx); + + break; // The connection face is found. Stop looping over the cell faces. Jump to next connection + } + } + } +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +const std::vector& RigNNCData::findConnectionIndices( size_t globalCellIndex, cvf::StructGridInterface::FaceType face) const +{ + ConnectionSearchMap::const_iterator it; + static std::vector empty; + + it = m_cellIdxToFaceToConnectionIdxMap.find(globalCellIndex); + if (it != m_cellIdxToFaceToConnectionIdxMap.end()) + { + return it->second[face]; + } + + return empty; } diff --git a/ApplicationCode/ReservoirDataModel/RigNNCData.h b/ApplicationCode/ReservoirDataModel/RigNNCData.h index cdefe31110..0b174ea426 100644 --- a/ApplicationCode/ReservoirDataModel/RigNNCData.h +++ b/ApplicationCode/ReservoirDataModel/RigNNCData.h @@ -35,22 +35,19 @@ class RigConnection public: RigConnection( ) : m_c1GlobIdx(cvf::UNDEFINED_SIZE_T), - m_c1Face(cvf::StructGridInterface::POS_I), + m_c1Face(cvf::StructGridInterface::NO_FACE), m_c2GlobIdx(cvf::UNDEFINED_SIZE_T), - m_c2Face(cvf::StructGridInterface::NEG_I), - m_hasNoSharedArea(false), m_transmissibility(0.0) {} - - - bool m_hasNoSharedArea; size_t m_c1GlobIdx; cvf::StructGridInterface::FaceType m_c1Face; size_t m_c2GlobIdx; - cvf::StructGridInterface::FaceType m_c2Face; //7 Probably Unused. Remove double m_transmissibility; + + std::vector m_polygon; + /* enum NNCType { @@ -64,12 +61,14 @@ class RigNNCData : public cvf::Object public: RigNNCData(); + const std::vector& findConnectionIndices(size_t globalCellIndex, cvf::StructGridInterface::FaceType face) const; void processConnections(const RigMainGrid& mainGrid); std::vector& connections() { return m_connections; } const std::vector& connections() const { return m_connections; }; private: - + typedef std::map, 7 > > ConnectionSearchMap; + ConnectionSearchMap m_cellIdxToFaceToConnectionIdxMap; std::vector m_connections; };