OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

(#540) WIP: Split the rangefilter mapping code into files

Browse Source 
Preparing for assembling it all.
This commit is contained in:
Jacob Støren 2015-10-19 11:27:56 +02:00
parent c023bc9b2a
commit a5ef83a2ee
7 changed files with 890 additions and 758 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
ApplicationCode/ReservoirDataModel/CMakeLists_filesNotToUnitTest.cmake
View File

@ -7,12 +7,16 @@ endif()
set (SOURCE_GROUP_HEADER_FILES set (SOURCE_GROUP_HEADER_FILES
${CEE_CURRENT_LIST_DIR}RigGeoMechWellLogExtractor.h ${CEE_CURRENT_LIST_DIR}RigGeoMechWellLogExtractor.h
${CEE_CURRENT_LIST_DIR}RigCaseToCaseCellMapper.h ${CEE_CURRENT_LIST_DIR}RigCaseToCaseCellMapper.h
${CEE_CURRENT_LIST_DIR}RigCaseToCaseCellMapperTools.h
${CEE_CURRENT_LIST_DIR}RigCaseToCaseRangeFilterMapper.h
${CEE_CURRENT_LIST_DIR}RigWellLogFile.h ${CEE_CURRENT_LIST_DIR}RigWellLogFile.h
) )
set (SOURCE_GROUP_SOURCE_FILES set (SOURCE_GROUP_SOURCE_FILES
${CEE_CURRENT_LIST_DIR}RigGeoMechWellLogExtractor.cpp ${CEE_CURRENT_LIST_DIR}RigGeoMechWellLogExtractor.cpp
${CEE_CURRENT_LIST_DIR}RigCaseToCaseCellMapper.cpp ${CEE_CURRENT_LIST_DIR}RigCaseToCaseCellMapper.cpp
${CEE_CURRENT_LIST_DIR}RigCaseToCaseCellMapperTools.cpp
${CEE_CURRENT_LIST_DIR}RigCaseToCaseRangeFilterMapper.cpp
${CEE_CURRENT_LIST_DIR}RigWellLogFile.cpp ${CEE_CURRENT_LIST_DIR}RigWellLogFile.cpp
) )

745
ApplicationCode/ReservoirDataModel/RigCaseToCaseCellMapper.cpp
View File

@ -18,6 +18,8 @@
///////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////
#include "RigCaseToCaseCellMapper.h" #include "RigCaseToCaseCellMapper.h"
#include "RigCaseToCaseCellMapperTools.h"
#include "RigFemPart.h" #include "RigFemPart.h"
#include "RigMainGrid.h" #include "RigMainGrid.h"
#include "RigFemPartGrid.h" #include "RigFemPartGrid.h"
@ -141,383 +143,6 @@ void RigCaseToCaseCellMapper::addMapping(int depCaseCellIdx, int masterCaseMatch
} }
//==================================================================================================
///
//==================================================================================================
class RigNeighborCornerFinder
{
public:
RigNeighborCornerFinder(const RigMainGrid* mainGrid, size_t baseI, size_t baseJ, size_t baseK)
: m_mainGrid(mainGrid),
m_baseI(baseI),
m_baseJ(baseJ),
m_baseK(baseK)
{}
const caf::SizeTArray8* neighborIndices(int offsetI, int offsetJ, int offsetK)
{
if (offsetI < 0 && m_baseI == 0) return NULL;
if (offsetJ < 0 && m_baseJ == 0) return NULL;
if (offsetK < 0 && m_baseK == 0) return NULL;
if (offsetI > 0 && m_baseI == m_mainGrid->cellCountI()-1) return NULL;
if (offsetJ > 0 && m_baseJ == m_mainGrid->cellCountJ()-1) return NULL;
if (offsetK > 0 && m_baseK == m_mainGrid->cellCountK()-1) return NULL;
size_t gridLocalCellIndex = m_mainGrid->cellIndexFromIJK(m_baseI + offsetI, m_baseJ + offsetJ, m_baseK + offsetK);
const RigCell& cell = m_mainGrid->cells()[gridLocalCellIndex];
return &(cell.cornerIndices());
}
private:
const RigMainGrid* m_mainGrid;
size_t m_baseI;
size_t m_baseJ;
size_t m_baseK;
};
//==================================================================================================
///
//==================================================================================================
//--------------------------------------------------------------------------------------------------
/// Average of neighbor corresponding nodes
//--------------------------------------------------------------------------------------------------
void RigCaseToCaseMapperTools::estimatedFemCellFromEclCell(const RigMainGrid* eclGrid, size_t reservoirCellIndex, cvf::Vec3d estimatedElmCorners[8])
{
CVF_TIGHT_ASSERT(reservoirCellIndex < eclGrid->cellCount()); // Assume reservoirCellIdx == localGridCellIdx for maingrid
const std::vector<cvf::Vec3d>& eclNodes = eclGrid->nodes();
size_t I,J,K;
eclGrid->ijkFromCellIndex(reservoirCellIndex, &I, &J, &K);
RigNeighborCornerFinder nbFinder(eclGrid, I,J,K);
// Cell corner Averaging mapping: Local cell index in neighbor matching specific corner of this cell
// N - Negative P - positive
// 0 <- NI[1] NINJ[2] NJ[3] NK[4] NINK[5] NINJNK[6] NJNK[7]
// 1 <- NJ[2] PINJ[3] PI[0] NK[5] NJNK[6] PINJNK[7] PINK[4]
// 2 <- PI[3] PIPJ[0] PJ[1] NK[6] PINK[7] PIPJNK[4] PJNK[5]
// 3 <- PJ[0] NIPJ[1] NI[2] NK[7] PJNK[4] NIPJNK[5] NINK[6]
// 4 <- NI[5] NINJ[6] NJ[7] PK[0] NIPK[1] NINJPK[2] NJPK[3]
// 5 <- NJ[6] PINJ[7] PI[4] PK[1] NJPK[2] PINJPK[3] PIPK[0]
// 6 <- PI[7] PIPJ[4] PJ[5] PK[2] PIPK[3] PIPJPK[0] PJPK[1]
// 7 <- PJ[4] NIPJ[5] NI[6] PK[3] PJPK[0] NIPJPK[1] NIPK[2]
const caf::SizeTArray8* IJK = nbFinder.neighborIndices( 0, 0, 0);
const caf::SizeTArray8* NI = nbFinder.neighborIndices(-1, 0, 0);
const caf::SizeTArray8* NJ = nbFinder.neighborIndices( 0,-1, 0);
const caf::SizeTArray8* PI = nbFinder.neighborIndices( 1, 0, 0);
const caf::SizeTArray8* PJ = nbFinder.neighborIndices( 0, 1, 0);
const caf::SizeTArray8* NK = nbFinder.neighborIndices( 0, 0,-1);
const caf::SizeTArray8* PK = nbFinder.neighborIndices( 0, 0, 1);
const caf::SizeTArray8* NINJ = nbFinder.neighborIndices(-1,-1, 0);
const caf::SizeTArray8* PINJ = nbFinder.neighborIndices( 1,-1, 0);
const caf::SizeTArray8* PIPJ = nbFinder.neighborIndices( 1, 1, 0);
const caf::SizeTArray8* NIPJ = nbFinder.neighborIndices(-1, 1, 0);
const caf::SizeTArray8* NINK = nbFinder.neighborIndices(-1, 0,-1);
const caf::SizeTArray8* NJNK = nbFinder.neighborIndices( 0,-1,-1);
const caf::SizeTArray8* PINK = nbFinder.neighborIndices( 1, 0,-1);
const caf::SizeTArray8* PJNK = nbFinder.neighborIndices( 0, 1,-1);
const caf::SizeTArray8* NIPK = nbFinder.neighborIndices(-1, 0, 1);
const caf::SizeTArray8* NJPK = nbFinder.neighborIndices( 0,-1, 1);
const caf::SizeTArray8* PIPK = nbFinder.neighborIndices( 1, 0, 1);
const caf::SizeTArray8* PJPK = nbFinder.neighborIndices( 0, 1, 1);
const caf::SizeTArray8* NINJNK = nbFinder.neighborIndices(-1,-1,-1);
const caf::SizeTArray8* PINJNK = nbFinder.neighborIndices( 1,-1,-1);
const caf::SizeTArray8* PIPJNK = nbFinder.neighborIndices( 1, 1,-1);
const caf::SizeTArray8* NIPJNK = nbFinder.neighborIndices(-1, 1,-1);
const caf::SizeTArray8* NINJPK = nbFinder.neighborIndices(-1,-1, 1);
const caf::SizeTArray8* PINJPK = nbFinder.neighborIndices( 1,-1, 1);
const caf::SizeTArray8* PIPJPK = nbFinder.neighborIndices( 1, 1, 1);
const caf::SizeTArray8* NIPJPK = nbFinder.neighborIndices(-1, 1, 1);
std::vector<size_t> contributingNodeIndicesPrCellCorner[8];
if (IJK ) contributingNodeIndicesPrCellCorner[0].push_back((*IJK )[0]);
if (NI ) contributingNodeIndicesPrCellCorner[0].push_back((*NI )[1]);
if (NINJ ) contributingNodeIndicesPrCellCorner[0].push_back((*NINJ )[2]);
if (NJ ) contributingNodeIndicesPrCellCorner[0].push_back((*NJ )[3]);
if (NK ) contributingNodeIndicesPrCellCorner[0].push_back((*NK )[4]);
if (NINK ) contributingNodeIndicesPrCellCorner[0].push_back((*NINK )[5]);
if (NINJNK) contributingNodeIndicesPrCellCorner[0].push_back((*NINJNK)[6]);
if (NJNK ) contributingNodeIndicesPrCellCorner[0].push_back((*NJNK )[7]);
if (IJK ) contributingNodeIndicesPrCellCorner[1].push_back((*IJK )[1]);
if (NJ ) contributingNodeIndicesPrCellCorner[1].push_back((*NJ )[2]);
if (PINJ ) contributingNodeIndicesPrCellCorner[1].push_back((*PINJ )[3]);
if (PI ) contributingNodeIndicesPrCellCorner[1].push_back((*PI )[0]);
if (NK ) contributingNodeIndicesPrCellCorner[1].push_back((*NK )[5]);
if (NJNK ) contributingNodeIndicesPrCellCorner[1].push_back((*NJNK )[6]);
if (PINJNK) contributingNodeIndicesPrCellCorner[1].push_back((*PINJNK)[7]);
if (PINK ) contributingNodeIndicesPrCellCorner[1].push_back((*PINK )[4]);
if (IJK ) contributingNodeIndicesPrCellCorner[2].push_back((*IJK )[2]);
if (PI ) contributingNodeIndicesPrCellCorner[2].push_back((*PI )[3]);
if (PIPJ ) contributingNodeIndicesPrCellCorner[2].push_back((*PIPJ )[0]);
if (PJ ) contributingNodeIndicesPrCellCorner[2].push_back((*PJ )[1]);
if (NK ) contributingNodeIndicesPrCellCorner[2].push_back((*NK )[6]);
if (PINK ) contributingNodeIndicesPrCellCorner[2].push_back((*PINK )[7]);
if (PIPJNK) contributingNodeIndicesPrCellCorner[2].push_back((*PIPJNK)[4]);
if (PJNK ) contributingNodeIndicesPrCellCorner[2].push_back((*PJNK )[5]);
if (IJK ) contributingNodeIndicesPrCellCorner[3].push_back((*IJK )[3]);
if (PJ ) contributingNodeIndicesPrCellCorner[3].push_back((*PJ )[0]);
if (NIPJ ) contributingNodeIndicesPrCellCorner[3].push_back((*NIPJ )[1]);
if (NI ) contributingNodeIndicesPrCellCorner[3].push_back((*NI )[2]);
if (NK ) contributingNodeIndicesPrCellCorner[3].push_back((*NK )[7]);
if (PJNK ) contributingNodeIndicesPrCellCorner[3].push_back((*PJNK )[4]);
if (NIPJNK) contributingNodeIndicesPrCellCorner[3].push_back((*NIPJNK)[5]);
if (NINK ) contributingNodeIndicesPrCellCorner[3].push_back((*NINK )[6]);
// 4 <- NI[5] NINJ[6] NJ[7] PK[0] NIPK[1] NINJPK[2] NJPK[3]
if (IJK ) contributingNodeIndicesPrCellCorner[4].push_back((*IJK )[4]);
if (NI ) contributingNodeIndicesPrCellCorner[4].push_back((*NI )[5]);
if (NINJ ) contributingNodeIndicesPrCellCorner[4].push_back((*NINJ )[6]);
if (NJ ) contributingNodeIndicesPrCellCorner[4].push_back((*NJ )[7]);
if (PK ) contributingNodeIndicesPrCellCorner[4].push_back((*PK )[0]);
if (NIPK ) contributingNodeIndicesPrCellCorner[4].push_back((*NIPK )[1]);
if (NINJPK) contributingNodeIndicesPrCellCorner[4].push_back((*NINJPK)[2]);
if (NJPK ) contributingNodeIndicesPrCellCorner[4].push_back((*NJPK )[3]);
if (IJK ) contributingNodeIndicesPrCellCorner[5].push_back((*IJK )[5]);
if (NJ ) contributingNodeIndicesPrCellCorner[5].push_back((*NJ )[6]);
if (PINJ ) contributingNodeIndicesPrCellCorner[5].push_back((*PINJ )[7]);
if (PI ) contributingNodeIndicesPrCellCorner[5].push_back((*PI )[4]);
if (PK ) contributingNodeIndicesPrCellCorner[5].push_back((*PK )[1]);
if (NJPK ) contributingNodeIndicesPrCellCorner[5].push_back((*NJPK )[2]);
if (PINJPK) contributingNodeIndicesPrCellCorner[5].push_back((*PINJPK)[3]);
if (PIPK ) contributingNodeIndicesPrCellCorner[5].push_back((*PIPK )[0]);
// 6 <- PI[7] PIPJ[4] PJ[5] PK[2] PIPK[3] PIPJPK[0] PJPK[1]
if (IJK ) contributingNodeIndicesPrCellCorner[6].push_back((*IJK )[6]);
if (PI ) contributingNodeIndicesPrCellCorner[6].push_back((*PI )[7]);
if (PIPJ ) contributingNodeIndicesPrCellCorner[6].push_back((*PIPJ )[4]);
if (PJ ) contributingNodeIndicesPrCellCorner[6].push_back((*PJ )[5]);
if (PK ) contributingNodeIndicesPrCellCorner[6].push_back((*PK )[2]);
if (PIPK ) contributingNodeIndicesPrCellCorner[6].push_back((*PIPK )[3]);
if (PIPJPK) contributingNodeIndicesPrCellCorner[6].push_back((*PIPJPK)[0]);
if (PJPK ) contributingNodeIndicesPrCellCorner[6].push_back((*PJPK )[1]);
if (IJK ) contributingNodeIndicesPrCellCorner[7].push_back((*IJK )[7]);
if (PJ ) contributingNodeIndicesPrCellCorner[7].push_back((*PJ )[4]);
if (NIPJ ) contributingNodeIndicesPrCellCorner[7].push_back((*NIPJ )[5]);
if (NI ) contributingNodeIndicesPrCellCorner[7].push_back((*NI )[6]);
if (PK ) contributingNodeIndicesPrCellCorner[7].push_back((*PK )[3]);
if (PJPK ) contributingNodeIndicesPrCellCorner[7].push_back((*PJPK )[0]);
if (NIPJPK) contributingNodeIndicesPrCellCorner[7].push_back((*NIPJPK)[1]);
if (NIPK ) contributingNodeIndicesPrCellCorner[7].push_back((*NIPK )[2]);
// Average the nodes
for (size_t cornIdx = 0; cornIdx < 8; ++cornIdx)
{
estimatedElmCorners[cornIdx] = cvf::Vec3d::ZERO;
size_t contribCount = contributingNodeIndicesPrCellCorner[cornIdx].size();
for (size_t ctnIdx = 0; ctnIdx < contribCount; ++ctnIdx)
{
estimatedElmCorners[cornIdx] += eclNodes[contributingNodeIndicesPrCellCorner[cornIdx][ctnIdx]];
}
estimatedElmCorners[cornIdx] /= contribCount;
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigCaseToCaseMapperTools::rotateQuad(cvf::Vec3d quad[4], int idxToNewStart)
{
if (idxToNewStart == 0) return;
cvf::Vec3d tmpQuad[4];
tmpQuad[0] = quad[0];
tmpQuad[1] = quad[1];
tmpQuad[2] = quad[2];
tmpQuad[3] = quad[3];
quad[0] = tmpQuad[idxToNewStart];
++idxToNewStart; if (idxToNewStart > 3) idxToNewStart = 0;
quad[1] = tmpQuad[idxToNewStart];
++idxToNewStart; if (idxToNewStart > 3) idxToNewStart = 0;
quad[2] = tmpQuad[idxToNewStart];
++idxToNewStart; if (idxToNewStart > 3) idxToNewStart = 0;
quad[3] = tmpQuad[idxToNewStart];
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigCaseToCaseMapperTools::flipQuadWinding(cvf::Vec3d quad[4])
{
cvf::Vec3d temp = quad[1];
quad[1] = quad[3];
quad[3] = temp;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
int RigCaseToCaseMapperTools::quadVxClosestToXYOfPoint( const cvf::Vec3d point, const cvf::Vec3d quad[4])
{
double minSqDist = HUGE_VAL;
int quadVxIdxClosestToPoint = cvf::UNDEFINED_INT;
for (int i = 0; i < 4; ++i)
{
cvf::Vec3d diff = quad[i]- point;
diff[2] = 0.0;
double sqDist = diff.lengthSquared();
if (sqDist < minSqDist)
{
minSqDist = sqDist;
quadVxIdxClosestToPoint = i;
}
}
return quadVxIdxClosestToPoint;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigCaseToCaseMapperTools::elementCorners(const RigFemPart* femPart, int elmIdx, cvf::Vec3d elmCorners[8])
{
RigElementType elmType = femPart->elementType(elmIdx);
if (!(elmType == HEX8 || elmType == HEX8P)) return false;
const std::vector<cvf::Vec3f>& nodeCoords = femPart->nodes().coordinates;
const int* cornerIndices = femPart->connectivities(elmIdx);
elmCorners[0] = cvf::Vec3d(nodeCoords[cornerIndices[0]]);
elmCorners[1] = cvf::Vec3d(nodeCoords[cornerIndices[1]]);
elmCorners[2] = cvf::Vec3d(nodeCoords[cornerIndices[2]]);
elmCorners[3] = cvf::Vec3d(nodeCoords[cornerIndices[3]]);
elmCorners[4] = cvf::Vec3d(nodeCoords[cornerIndices[4]]);
elmCorners[5] = cvf::Vec3d(nodeCoords[cornerIndices[5]]);
elmCorners[6] = cvf::Vec3d(nodeCoords[cornerIndices[6]]);
elmCorners[7] = cvf::Vec3d(nodeCoords[cornerIndices[7]]);
return true;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
int RigCaseToCaseMapperTools::findMatchingPOSKFaceIdx(const cvf::Vec3d baseCell[8],bool isBaseCellNormalsOutwards, const cvf::Vec3d c2[8])
{
int faceNodeCount;
const int* posKFace = RigFemTypes::localElmNodeIndicesForFace(HEX8, (int)(cvf::StructGridInterface::POS_K), &faceNodeCount);
double sign = isBaseCellNormalsOutwards ? 1.0 : -1.0;
cvf::Vec3d posKnormal = sign*(baseCell[posKFace[2]] - baseCell[posKFace[0]]) ^ (baseCell[posKFace[3]] - baseCell[posKFace[1]]);
posKnormal.normalize();
double minDiff = HUGE_VAL;
int bestFace = -1;
for (int faceIdx = 5; faceIdx >= 0; --faceIdx) // Backwards. might hit earlier more often
{
const int* face = RigFemTypes::localElmNodeIndicesForFace(HEX8, faceIdx, &faceNodeCount);
cvf::Vec3d normal = (c2[face[2]] - c2[face[0]]) ^ (c2[face[3]] - c2[face[1]]);
normal.normalize();
double sqDiff = (posKnormal-normal).lengthSquared();
if (sqDiff < minDiff)
{
minDiff = sqDiff;
bestFace = faceIdx;
if (minDiff < 0.1*0.1) break; // This must be the one. Do not search further
}
}
return bestFace;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigCaseToCaseMapperTools::isEclFemCellsMatching(const cvf::Vec3d baseCell[8],
cvf::Vec3d cell[8],
double xyTolerance, double zTolerance)
{
bool isMatching = true;
for (int i = 0; i < 4 ; ++i)
{
cvf::Vec3d diff = cell[i] - baseCell[i];
if (!(fabs(diff.x()) < xyTolerance && fabs(diff.y()) < xyTolerance && fabs(diff.z()) < zTolerance))
{
isMatching = false;
break;
}
}
return isMatching;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigCaseToCaseMapperTools::rotateCellTopologicallyToMatchBaseCell(const cvf::Vec3d * baseCell, bool baseCellFaceNormalsIsOutwards, cvf::Vec3d * cell)
{
int femDeepZFaceIdx = findMatchingPOSKFaceIdx(baseCell, baseCellFaceNormalsIsOutwards, cell);
{
cvf::Vec3d tmpFemCorners[8];
tmpFemCorners[0] = cell[0];
tmpFemCorners[1] = cell[1];
tmpFemCorners[2] = cell[2];
tmpFemCorners[3] = cell[3];
tmpFemCorners[4] = cell[4];
tmpFemCorners[5] = cell[5];
tmpFemCorners[6] = cell[6];
tmpFemCorners[7] = cell[7];
int femShallowZFaceIdx = RigFemTypes::oppositeFace(HEX8, femDeepZFaceIdx);
int faceNodeCount;
const int* localElmNodeIndicesForPOSKFace = RigFemTypes::localElmNodeIndicesForFace(HEX8, femDeepZFaceIdx, &faceNodeCount);
const int* localElmNodeIndicesForNEGKFace = RigFemTypes::localElmNodeIndicesForFace(HEX8, femShallowZFaceIdx, &faceNodeCount);
cell[0] = tmpFemCorners[localElmNodeIndicesForNEGKFace[0]];
cell[1] = tmpFemCorners[localElmNodeIndicesForNEGKFace[1]];
cell[2] = tmpFemCorners[localElmNodeIndicesForNEGKFace[2]];
cell[3] = tmpFemCorners[localElmNodeIndicesForNEGKFace[3]];
cell[4] = tmpFemCorners[localElmNodeIndicesForPOSKFace[0]];
cell[5] = tmpFemCorners[localElmNodeIndicesForPOSKFace[1]];
cell[6] = tmpFemCorners[localElmNodeIndicesForPOSKFace[2]];
cell[7] = tmpFemCorners[localElmNodeIndicesForPOSKFace[3]];
}
cvf::Vec3d* femDeepestQuad = &(cell[4]);
cvf::Vec3d* femShallowQuad = &(cell[0]);
// Now the top/bottom have opposite winding. To make the comparisons and index rotations simpler
// flip the winding of the top or bottom face depending on whether the eclipse grid is inside-out
if (baseCellFaceNormalsIsOutwards)
{
flipQuadWinding(femShallowQuad);
}
else
{
flipQuadWinding(femDeepestQuad);
}
// We now need to rotate the fem quads to be alligned with the ecl quads
// Since the start point of the quad always is aligned with the opposite face-quad start
// we can find the rotation for the top, and apply it to both top and bottom
int femQuadStartIdx = quadVxClosestToXYOfPoint(baseCell[0], femShallowQuad);
rotateQuad(femDeepestQuad, femQuadStartIdx);
rotateQuad(femShallowQuad, femQuadStartIdx);
}
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
/// ///
//-------------------------------------------------------------------------------------------------- //--------------------------------------------------------------------------------------------------
@ -548,7 +173,7 @@ void RigCaseToCaseCellMapper::calculateEclToGeomCellMapping(RigMainGrid* masterE
#endif #endif
cvf::Vec3d geoMechConvertedEclCell[8]; cvf::Vec3d geoMechConvertedEclCell[8];
RigCaseToCaseMapperTools::estimatedFemCellFromEclCell(masterEclGrid, cellIdx, geoMechConvertedEclCell); RigCaseToCaseCellMapperTools::estimatedFemCellFromEclCell(masterEclGrid, cellIdx, geoMechConvertedEclCell);
cvf::BoundingBox elmBBox; cvf::BoundingBox elmBBox;
for (int i = 0; i < 8 ; ++i) elmBBox.add(geoMechConvertedEclCell[i]); for (int i = 0; i < 8 ; ++i) elmBBox.add(geoMechConvertedEclCell[i]);
@ -560,11 +185,11 @@ void RigCaseToCaseCellMapper::calculateEclToGeomCellMapping(RigMainGrid* masterE
{ {
int elmIdx = static_cast<int>(closeElements[ccIdx]); int elmIdx = static_cast<int>(closeElements[ccIdx]);
RigCaseToCaseMapperTools::elementCorners(dependentFemPart, elmIdx, elmCorners); RigCaseToCaseCellMapperTools::elementCorners(dependentFemPart, elmIdx, elmCorners);
RigCaseToCaseMapperTools::rotateCellTopologicallyToMatchBaseCell(geoMechConvertedEclCell, isEclFaceNormalsOutwards , elmCorners); RigCaseToCaseCellMapperTools::rotateCellTopologicallyToMatchBaseCell(geoMechConvertedEclCell, isEclFaceNormalsOutwards , elmCorners);
bool isMatching = RigCaseToCaseMapperTools::isEclFemCellsMatching(geoMechConvertedEclCell, elmCorners, bool isMatching = RigCaseToCaseCellMapperTools::isEclFemCellsMatching(geoMechConvertedEclCell, elmCorners,
xyTolerance, zTolerance); xyTolerance, zTolerance);
if (isMatching) if (isMatching)
@ -579,361 +204,3 @@ void RigCaseToCaseCellMapper::calculateEclToGeomCellMapping(RigMainGrid* masterE
} }
} }
} }
#if 0 // Inside Bounding box test
cvf::BoundingBox cellBBox;
for (int i = 0; i < 8 ; ++i) cellBBox.add(cellCorners[i]);
cvf::Vec3d cs = cellBBox.min();
cvf::Vec3d cl = cellBBox.max();
cvf::Vec3d es = elmBBox.min();
cvf::Vec3d el = elmBBox.max();
if ( ( (cs.x() + xyTolerance) >= es.x() && (cl.x() - xyTolerance) <= el.x())
&& ( (cs.y() + xyTolerance) >= es.y() && (cl.y() - xyTolerance) <= el.y())
&& ( (cs.z() + zTolerance ) >= es.z() && (cl.z() - zTolerance ) <= el.z()) )
{
// Cell bb equal or inside Elm bb
isMatching = true;
}
if ( ( (es.x() + xyTolerance) >= cs.x() && (el.x() - xyTolerance) <= cl.x())
&& ( (es.y() + xyTolerance) >= cs.y() && (el.y() - xyTolerance) <= cl.y())
&& ( (es.z() + zTolerance ) >= cs.z() && (el.z() - zTolerance ) <= cl.z()) )
{
// Elm bb equal or inside Cell bb
isMatching = true;
}
#endif
#if 0
{
const std::vector<cvf::Vec3d>& eclNodes = eclGrid->nodes();
const RigCell& cell = eclGrid->cells()[reservoirCellIndex];
const caf::SizeTArray8& cornerIndices = cell.cornerIndices();
int faceNodeCount;
const int* localElmNodeIndicesForTopZFace = RigFemTypes::localElmNodeIndicesForFace(HEX8, 4, &faceNodeCount);
const int* localElmNodeIndicesForBotZFace = RigFemTypes::localElmNodeIndicesForFace(HEX8, 5, &faceNodeCount);
eclDeepestQuad[0] = eclNodes[cornerIndices[localElmNodeIndicesForTopZFace[0]]];
eclDeepestQuad[1] = eclNodes[cornerIndices[localElmNodeIndicesForTopZFace[1]]];
eclDeepestQuad[2] = eclNodes[cornerIndices[localElmNodeIndicesForTopZFace[2]]];
eclDeepestQuad[3] = eclNodes[cornerIndices[localElmNodeIndicesForTopZFace[3]]];
eclShallowQuad[0] = eclNodes[cornerIndices[localElmNodeIndicesForBotZFace[0]]];
eclShallowQuad[1] = eclNodes[cornerIndices[localElmNodeIndicesForBotZFace[1]]];
eclShallowQuad[2] = eclNodes[cornerIndices[localElmNodeIndicesForBotZFace[2]]];
eclShallowQuad[3] = eclNodes[cornerIndices[localElmNodeIndicesForBotZFace[3]]];
}
#endif
#if 0
// First search K=1 diagonally for a seed cell; A cell without collapsings, and without faults
size_t minIJCount = masterEclGrid->cellCountI();
if (minIJCount > masterEclGrid->cellCountJ())
minIJCount = masterEclGrid->cellCountJ();
for (size_t ij = 0; ij < minIJCount; ++ij )
{
size_t localCellIdx = masterEclGrid->cellIndexFromIJK(ij, ij, 0);
size_t reservoirCellIdx = masterEclGrid->reservoirCellIndex(localCellIdx);
cvf::Vec3d vertices[8];
masterEclGrid->cellCornerVertices(localCellIdx, vertices);
if (!isCellNormal(vertices))
continue;
const RigFault* fault = masterEclGrid->findFaultFromCellIndexAndCellFace(reservoirCellIdx, cvf::StructGridInterface::POS_I);
}
#endif
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigCaseToCaseMapperTools::calculateCellCenter(cvf::Vec3d elmCorners[8])
{
cvf::Vec3d avg(cvf::Vec3d::ZERO);
size_t i;
for (i = 0; i < 8; i++)
{
avg += elmCorners[i];
}
avg /= 8.0;
return avg;
}
#include "RimCellRangeFilter.h"
class RigCaseToCaseRangeFilterMapper
{
public:
static void convertRangeFilterEclToFem(RimCellRangeFilter* srcFilter, const RigMainGrid* srcEclGrid,
RimCellRangeFilter* dstFilter, const RigFemPart* dstFemPart)
{
CVF_ASSERT(srcFilter && srcEclGrid && dstFilter && dstFemPart);
CVF_ASSERT(srcFilter->gridIndex() == 0); // LGR not supported yet
struct RangeFilterCorner { RangeFilterCorner(): isExactMatch(false){} cvf::Vec3st ijk; bool isExactMatch; };
RangeFilterCorner rangeFilterMatches[8];
size_t srcStartI = srcFilter->startIndexI() - 1;
size_t srcStartJ = srcFilter->startIndexJ() - 1;
size_t srcStartK = srcFilter->startIndexK() - 1;
size_t srcEndI = srcStartI + srcFilter->cellCountI();
size_t srcEndJ = srcStartJ + srcFilter->cellCountJ();
size_t srcEndK = srcStartK + srcFilter->cellCountK();
cvf::Vec3st srcRangeCube[8];
srcRangeCube[0] = cvf::Vec3st(srcStartI, srcStartJ, srcStartK);
srcRangeCube[1] = cvf::Vec3st(srcEndI, srcStartJ, srcStartK);
srcRangeCube[2] = cvf::Vec3st(srcEndI, srcEndJ, srcStartK);
srcRangeCube[3] = cvf::Vec3st(srcStartI, srcEndJ, srcStartK);
srcRangeCube[4] = cvf::Vec3st(srcStartI, srcStartJ, srcEndK);
srcRangeCube[5] = cvf::Vec3st(srcEndI, srcStartJ, srcEndK);
srcRangeCube[6] = cvf::Vec3st(srcEndI, srcEndJ, srcEndK);
srcRangeCube[7] = cvf::Vec3st(srcStartI, srcEndJ, srcEndK);
size_t dstStartI = cvf::UNDEFINED_SIZE_T;
size_t dstStartJ = cvf::UNDEFINED_SIZE_T;
size_t dstStartK = cvf::UNDEFINED_SIZE_T;
size_t dstEndI = cvf::UNDEFINED_SIZE_T;
size_t dstEndJ = cvf::UNDEFINED_SIZE_T;
size_t dstEndK = cvf::UNDEFINED_SIZE_T;
bool foundExactMatch = false;
int cornerIdx = 0;
int diagIdx = 6;// Index to diagonal corner
for (cornerIdx = 0; cornerIdx < 4; ++cornerIdx)
{
diagIdx = (cornerIdx < 2) ? cornerIdx + 6 : cornerIdx + 2;
rangeFilterMatches[cornerIdx].isExactMatch = findBestFemCellFromEclCell(srcEclGrid,
srcRangeCube[cornerIdx][0],
srcRangeCube[cornerIdx][1],
srcRangeCube[cornerIdx][2],
dstFemPart,
&(rangeFilterMatches[cornerIdx].ijk[0]),
&(rangeFilterMatches[cornerIdx].ijk[1]),
&(rangeFilterMatches[cornerIdx].ijk[2]));
rangeFilterMatches[diagIdx].isExactMatch = findBestFemCellFromEclCell(srcEclGrid,
srcRangeCube[diagIdx][0],
srcRangeCube[diagIdx][1],
srcRangeCube[diagIdx][2],
dstFemPart,
&(rangeFilterMatches[diagIdx].ijk[0]),
&(rangeFilterMatches[diagIdx].ijk[1]),
&(rangeFilterMatches[diagIdx].ijk[2]));
if (rangeFilterMatches[cornerIdx].isExactMatch && rangeFilterMatches[diagIdx].isExactMatch)
{
foundExactMatch = true;
break;
}
}
// Get the start and end IJK from the matched corners
if (foundExactMatch)
{
// Populate dst range filter from the diagonal that matches exact
dstStartI = CVF_MIN(rangeFilterMatches[cornerIdx].ijk[0], rangeFilterMatches[diagIdx].ijk[0]);
dstStartJ = CVF_MIN(rangeFilterMatches[cornerIdx].ijk[1], rangeFilterMatches[diagIdx].ijk[1]);
dstStartK = CVF_MIN(rangeFilterMatches[cornerIdx].ijk[2], rangeFilterMatches[diagIdx].ijk[2]);
dstEndI = CVF_MAX(rangeFilterMatches[cornerIdx].ijk[0], rangeFilterMatches[diagIdx].ijk[0]);
dstEndJ = CVF_MAX(rangeFilterMatches[cornerIdx].ijk[1], rangeFilterMatches[diagIdx].ijk[1]);
dstEndK = CVF_MAX(rangeFilterMatches[cornerIdx].ijk[2], rangeFilterMatches[diagIdx].ijk[2]);
}
else
{
// Todo: be even smarter, and use possible matching corners to add up an as best solution as possible.
// For now we just take the first diagonal.
dstStartI = rangeFilterMatches[0].ijk[0];
dstStartJ = rangeFilterMatches[0].ijk[1];
dstStartK = rangeFilterMatches[0].ijk[2];
dstEndI = rangeFilterMatches[6].ijk[0];
dstEndJ = rangeFilterMatches[6].ijk[1];
dstEndK = rangeFilterMatches[6].ijk[2];
}
// Populate the dst range filter with new data
if ( (dstStartI != cvf::UNDEFINED_SIZE_T && dstStartJ != cvf::UNDEFINED_SIZE_T && dstStartK != cvf::UNDEFINED_SIZE_T)
&& (dstEndI != cvf::UNDEFINED_SIZE_T && dstEndJ != cvf::UNDEFINED_SIZE_T && dstEndK != cvf::UNDEFINED_SIZE_T))
{
dstFilter->startIndexJ = static_cast<int>(dstStartI + 1);
dstFilter->startIndexK = static_cast<int>(dstStartJ + 1);
dstFilter->startIndexI = static_cast<int>(dstStartK + 1);
dstFilter->cellCountI = static_cast<int>(dstEndI - (dstStartI-1));
dstFilter->cellCountJ = static_cast<int>(dstEndJ - (dstStartJ-1));
dstFilter->cellCountK = static_cast<int>(dstEndK - (dstStartK-1));
}
else
{
dstFilter->startIndexI = 1;
dstFilter->startIndexJ = 1;
dstFilter->startIndexK = 1;
dstFilter->cellCountI = 0;
dstFilter->cellCountJ = 0;
dstFilter->cellCountK = 0;
dstFilter->computeAndSetValidValues();
}
}
static bool findBestFemCellFromEclCell(const RigMainGrid* masterEclGrid, size_t ei, size_t ej, size_t ek,
const RigFemPart* dependentFemPart, size_t* fi, size_t * fj, size_t* fk)
{
// Find tolerance
double cellSizeI, cellSizeJ, cellSizeK;
masterEclGrid->characteristicCellSizes(&cellSizeI, &cellSizeJ, &cellSizeK);
double xyTolerance = cellSizeI* 0.4;
double zTolerance = cellSizeK* 0.4;
bool isEclFaceNormalsOutwards = masterEclGrid->isFaceNormalsOutwards();
size_t cellIdx = masterEclGrid->cellIndexFromIJK( ei, ej, ek);
cvf::Vec3d geoMechConvertedEclCell[8];
RigCaseToCaseMapperTools::estimatedFemCellFromEclCell(masterEclGrid, cellIdx, geoMechConvertedEclCell);
cvf::BoundingBox elmBBox;
for (int i = 0; i < 8 ; ++i) elmBBox.add(geoMechConvertedEclCell[i]);
std::vector<size_t> closeElements;
dependentFemPart->findIntersectingCells(elmBBox, &closeElements);
cvf::Vec3d elmCorners[8];
int elmIdxToBestMatch = -1;
double sqDistToClosestElmCenter = HUGE_VAL;
cvf::Vec3d convEclCellCenter = RigCaseToCaseMapperTools::calculateCellCenter(geoMechConvertedEclCell);
bool foundExactMatch = false;
for (size_t ccIdx = 0; ccIdx < closeElements.size(); ++ccIdx)
{
int elmIdx = static_cast<int>(closeElements[ccIdx]);
RigCaseToCaseMapperTools::elementCorners(dependentFemPart, elmIdx, elmCorners);
cvf::Vec3d cellCenter = RigCaseToCaseMapperTools::calculateCellCenter(elmCorners);
double sqDist = (cellCenter - convEclCellCenter).lengthSquared();
if (sqDist < sqDistToClosestElmCenter)
{
elmIdxToBestMatch = elmIdx;
sqDistToClosestElmCenter = sqDist;
}
RigCaseToCaseMapperTools::rotateCellTopologicallyToMatchBaseCell(geoMechConvertedEclCell, isEclFaceNormalsOutwards, elmCorners);
foundExactMatch = RigCaseToCaseMapperTools::isEclFemCellsMatching(geoMechConvertedEclCell, elmCorners,
xyTolerance, zTolerance);
if (foundExactMatch)
{
elmIdxToBestMatch = elmIdx;
break;
}
}
if (elmIdxToBestMatch != -1)
{
dependentFemPart->structGrid()->ijkFromCellIndex(elmIdxToBestMatch, fi, fj, fk);
}
else
{
(*fi) = cvf::UNDEFINED_SIZE_T;
(*fj) = cvf::UNDEFINED_SIZE_T;
(*fk) = cvf::UNDEFINED_SIZE_T;
}
return foundExactMatch;
}
static bool findBestEclCellFromFemCell(const RigFemPart* dependentFemPart, size_t fi, size_t fj, size_t fk,
const RigMainGrid* masterEclGrid, size_t* ei, size_t* ej, size_t* ek)
{
// Find tolerance
double cellSizeI, cellSizeJ, cellSizeK;
masterEclGrid->characteristicCellSizes(&cellSizeI, &cellSizeJ, &cellSizeK);
double xyTolerance = cellSizeI* 0.4;
double zTolerance = cellSizeK* 0.4;
bool isEclFaceNormalsOutwards = masterEclGrid->isFaceNormalsOutwards();
int elementIdx = static_cast<int>(dependentFemPart->structGrid()->cellIndexFromIJK( fi, fj, fk));
cvf::Vec3d elmCorners[8];
RigCaseToCaseMapperTools::elementCorners(dependentFemPart, elementIdx, elmCorners);
cvf::BoundingBox elmBBox;
for (int i = 0; i < 8 ; ++i) elmBBox.add(elmCorners[i]);
std::vector<size_t> closeCells;
masterEclGrid->findIntersectingCells(elmBBox, &closeCells); // This might actually miss the exact one, but we have no other alternative yet.
size_t globCellIdxToBestMatch = cvf::UNDEFINED_SIZE_T;
double sqDistToClosestCellCenter = HUGE_VAL;
cvf::Vec3d elmCenter = RigCaseToCaseMapperTools::calculateCellCenter(elmCorners);
bool foundExactMatch = false;
for (size_t ccIdx = 0; ccIdx < closeCells.size(); ++ccIdx)
{
size_t cellIdx = closeCells[ccIdx];
cvf::Vec3d geoMechConvertedEclCell[8];
RigCaseToCaseMapperTools::estimatedFemCellFromEclCell(masterEclGrid, cellIdx, geoMechConvertedEclCell);
cvf::Vec3d cellCenter = RigCaseToCaseMapperTools::calculateCellCenter(geoMechConvertedEclCell);
double sqDist = (cellCenter - elmCenter).lengthSquared();
if (sqDist < sqDistToClosestCellCenter)
{
globCellIdxToBestMatch = cellIdx;
sqDistToClosestCellCenter = sqDist;
}
RigCaseToCaseMapperTools::rotateCellTopologicallyToMatchBaseCell(geoMechConvertedEclCell, isEclFaceNormalsOutwards, elmCorners);
foundExactMatch = RigCaseToCaseMapperTools::isEclFemCellsMatching(geoMechConvertedEclCell, elmCorners,
xyTolerance, zTolerance);
if (foundExactMatch)
{
globCellIdxToBestMatch = cellIdx;
break;
}
}
if (globCellIdxToBestMatch != cvf::UNDEFINED_SIZE_T)
{
masterEclGrid->ijkFromCellIndex(globCellIdxToBestMatch, ei, ej, ek);
}
else
{
(*ei) = cvf::UNDEFINED_SIZE_T;
(*ej) = cvf::UNDEFINED_SIZE_T;
(*ek) = cvf::UNDEFINED_SIZE_T;
}
return foundExactMatch;
}
};

19
ApplicationCode/ReservoirDataModel/RigCaseToCaseCellMapper.h
View File

@ -63,22 +63,3 @@ private:
RigFemPart* m_dependentFemPart; RigFemPart* m_dependentFemPart;
}; };
//==================================================================================================
///
//==================================================================================================
class RigCaseToCaseMapperTools
{
public:
static void estimatedFemCellFromEclCell(const RigMainGrid* eclGrid, size_t reservoirCellIndex, cvf::Vec3d estimatedElmCorners[8]);
static void rotateQuad(cvf::Vec3d quad[4], int idxToNewStart);
static void flipQuadWinding(cvf::Vec3d quad[4]);
static int quadVxClosestToXYOfPoint(const cvf::Vec3d point, const cvf::Vec3d quad[4]);
static bool elementCorners(const RigFemPart* femPart, int elmIdx, cvf::Vec3d elmCorners[8]);
static int findMatchingPOSKFaceIdx(const cvf::Vec3d baseCell[8], bool isBaseCellNormalsOutwards, const cvf::Vec3d c2[8]);
static bool isEclFemCellsMatching(const cvf::Vec3d baseCell[8], cvf::Vec3d cell[8], double xyTolerance, double zTolerance);
static void rotateCellTopologicallyToMatchBaseCell(const cvf::Vec3d * baseCell, bool baseCellFaceNormalsIsOutwards, cvf::Vec3d * cell);
static cvf::Vec3d calculateCellCenter(cvf::Vec3d elmCorners[8]);
};

492
ApplicationCode/ReservoirDataModel/RigCaseToCaseCellMapperTools.cpp Normal file
View File

@ -0,0 +1,492 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2015- Statoil ASA
// Copyright (C) 2015- Ceetron Solutions AS
//
// 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 "RigCaseToCaseCellMapper.h"
#include "RigCaseToCaseCellMapperTools.h"
#include "RigFemPart.h"
#include "RigMainGrid.h"
#include "RigFemPartGrid.h"
//==================================================================================================
///
//==================================================================================================
class RigNeighborCornerFinder
{
public:
RigNeighborCornerFinder(const RigMainGrid* mainGrid, size_t baseI, size_t baseJ, size_t baseK)
: m_mainGrid(mainGrid),
m_baseI(baseI),
m_baseJ(baseJ),
m_baseK(baseK)
{}
const caf::SizeTArray8* neighborIndices(int offsetI, int offsetJ, int offsetK)
{
if (offsetI < 0 && m_baseI == 0) return NULL;
if (offsetJ < 0 && m_baseJ == 0) return NULL;
if (offsetK < 0 && m_baseK == 0) return NULL;
if (offsetI > 0 && m_baseI == m_mainGrid->cellCountI()-1) return NULL;
if (offsetJ > 0 && m_baseJ == m_mainGrid->cellCountJ()-1) return NULL;
if (offsetK > 0 && m_baseK == m_mainGrid->cellCountK()-1) return NULL;
size_t gridLocalCellIndex = m_mainGrid->cellIndexFromIJK(m_baseI + offsetI, m_baseJ + offsetJ, m_baseK + offsetK);
const RigCell& cell = m_mainGrid->cells()[gridLocalCellIndex];
return &(cell.cornerIndices());
}
private:
const RigMainGrid* m_mainGrid;
size_t m_baseI;
size_t m_baseJ;
size_t m_baseK;
};
//==================================================================================================
///
//==================================================================================================
//--------------------------------------------------------------------------------------------------
/// Average of neighbor corresponding nodes
//--------------------------------------------------------------------------------------------------
void RigCaseToCaseCellMapperTools::estimatedFemCellFromEclCell(const RigMainGrid* eclGrid, size_t reservoirCellIndex, cvf::Vec3d estimatedElmCorners[8])
{
CVF_TIGHT_ASSERT(reservoirCellIndex < eclGrid->cellCount()); // Assume reservoirCellIdx == localGridCellIdx for maingrid
const std::vector<cvf::Vec3d>& eclNodes = eclGrid->nodes();
size_t I,J,K;
eclGrid->ijkFromCellIndex(reservoirCellIndex, &I, &J, &K);
RigNeighborCornerFinder nbFinder(eclGrid, I,J,K);
// Cell corner Averaging mapping: Local cell index in neighbor matching specific corner of this cell
// N - Negative P - positive
// 0 <- NI[1] NINJ[2] NJ[3] NK[4] NINK[5] NINJNK[6] NJNK[7]
// 1 <- NJ[2] PINJ[3] PI[0] NK[5] NJNK[6] PINJNK[7] PINK[4]
// 2 <- PI[3] PIPJ[0] PJ[1] NK[6] PINK[7] PIPJNK[4] PJNK[5]
// 3 <- PJ[0] NIPJ[1] NI[2] NK[7] PJNK[4] NIPJNK[5] NINK[6]
// 4 <- NI[5] NINJ[6] NJ[7] PK[0] NIPK[1] NINJPK[2] NJPK[3]
// 5 <- NJ[6] PINJ[7] PI[4] PK[1] NJPK[2] PINJPK[3] PIPK[0]
// 6 <- PI[7] PIPJ[4] PJ[5] PK[2] PIPK[3] PIPJPK[0] PJPK[1]
// 7 <- PJ[4] NIPJ[5] NI[6] PK[3] PJPK[0] NIPJPK[1] NIPK[2]
const caf::SizeTArray8* IJK = nbFinder.neighborIndices( 0, 0, 0);
const caf::SizeTArray8* NI = nbFinder.neighborIndices(-1, 0, 0);
const caf::SizeTArray8* NJ = nbFinder.neighborIndices( 0,-1, 0);
const caf::SizeTArray8* PI = nbFinder.neighborIndices( 1, 0, 0);
const caf::SizeTArray8* PJ = nbFinder.neighborIndices( 0, 1, 0);
const caf::SizeTArray8* NK = nbFinder.neighborIndices( 0, 0,-1);
const caf::SizeTArray8* PK = nbFinder.neighborIndices( 0, 0, 1);
const caf::SizeTArray8* NINJ = nbFinder.neighborIndices(-1,-1, 0);
const caf::SizeTArray8* PINJ = nbFinder.neighborIndices( 1,-1, 0);
const caf::SizeTArray8* PIPJ = nbFinder.neighborIndices( 1, 1, 0);
const caf::SizeTArray8* NIPJ = nbFinder.neighborIndices(-1, 1, 0);
const caf::SizeTArray8* NINK = nbFinder.neighborIndices(-1, 0,-1);
const caf::SizeTArray8* NJNK = nbFinder.neighborIndices( 0,-1,-1);
const caf::SizeTArray8* PINK = nbFinder.neighborIndices( 1, 0,-1);
const caf::SizeTArray8* PJNK = nbFinder.neighborIndices( 0, 1,-1);
const caf::SizeTArray8* NIPK = nbFinder.neighborIndices(-1, 0, 1);
const caf::SizeTArray8* NJPK = nbFinder.neighborIndices( 0,-1, 1);
const caf::SizeTArray8* PIPK = nbFinder.neighborIndices( 1, 0, 1);
const caf::SizeTArray8* PJPK = nbFinder.neighborIndices( 0, 1, 1);
const caf::SizeTArray8* NINJNK = nbFinder.neighborIndices(-1,-1,-1);
const caf::SizeTArray8* PINJNK = nbFinder.neighborIndices( 1,-1,-1);
const caf::SizeTArray8* PIPJNK = nbFinder.neighborIndices( 1, 1,-1);
const caf::SizeTArray8* NIPJNK = nbFinder.neighborIndices(-1, 1,-1);
const caf::SizeTArray8* NINJPK = nbFinder.neighborIndices(-1,-1, 1);
const caf::SizeTArray8* PINJPK = nbFinder.neighborIndices( 1,-1, 1);
const caf::SizeTArray8* PIPJPK = nbFinder.neighborIndices( 1, 1, 1);
const caf::SizeTArray8* NIPJPK = nbFinder.neighborIndices(-1, 1, 1);
std::vector<size_t> contributingNodeIndicesPrCellCorner[8];
if (IJK ) contributingNodeIndicesPrCellCorner[0].push_back((*IJK )[0]);
if (NI ) contributingNodeIndicesPrCellCorner[0].push_back((*NI )[1]);
if (NINJ ) contributingNodeIndicesPrCellCorner[0].push_back((*NINJ )[2]);
if (NJ ) contributingNodeIndicesPrCellCorner[0].push_back((*NJ )[3]);
if (NK ) contributingNodeIndicesPrCellCorner[0].push_back((*NK )[4]);
if (NINK ) contributingNodeIndicesPrCellCorner[0].push_back((*NINK )[5]);
if (NINJNK) contributingNodeIndicesPrCellCorner[0].push_back((*NINJNK)[6]);
if (NJNK ) contributingNodeIndicesPrCellCorner[0].push_back((*NJNK )[7]);
if (IJK ) contributingNodeIndicesPrCellCorner[1].push_back((*IJK )[1]);
if (NJ ) contributingNodeIndicesPrCellCorner[1].push_back((*NJ )[2]);
if (PINJ ) contributingNodeIndicesPrCellCorner[1].push_back((*PINJ )[3]);
if (PI ) contributingNodeIndicesPrCellCorner[1].push_back((*PI )[0]);
if (NK ) contributingNodeIndicesPrCellCorner[1].push_back((*NK )[5]);
if (NJNK ) contributingNodeIndicesPrCellCorner[1].push_back((*NJNK )[6]);
if (PINJNK) contributingNodeIndicesPrCellCorner[1].push_back((*PINJNK)[7]);
if (PINK ) contributingNodeIndicesPrCellCorner[1].push_back((*PINK )[4]);
if (IJK ) contributingNodeIndicesPrCellCorner[2].push_back((*IJK )[2]);
if (PI ) contributingNodeIndicesPrCellCorner[2].push_back((*PI )[3]);
if (PIPJ ) contributingNodeIndicesPrCellCorner[2].push_back((*PIPJ )[0]);
if (PJ ) contributingNodeIndicesPrCellCorner[2].push_back((*PJ )[1]);
if (NK ) contributingNodeIndicesPrCellCorner[2].push_back((*NK )[6]);
if (PINK ) contributingNodeIndicesPrCellCorner[2].push_back((*PINK )[7]);
if (PIPJNK) contributingNodeIndicesPrCellCorner[2].push_back((*PIPJNK)[4]);
if (PJNK ) contributingNodeIndicesPrCellCorner[2].push_back((*PJNK )[5]);
if (IJK ) contributingNodeIndicesPrCellCorner[3].push_back((*IJK )[3]);
if (PJ ) contributingNodeIndicesPrCellCorner[3].push_back((*PJ )[0]);
if (NIPJ ) contributingNodeIndicesPrCellCorner[3].push_back((*NIPJ )[1]);
if (NI ) contributingNodeIndicesPrCellCorner[3].push_back((*NI )[2]);
if (NK ) contributingNodeIndicesPrCellCorner[3].push_back((*NK )[7]);
if (PJNK ) contributingNodeIndicesPrCellCorner[3].push_back((*PJNK )[4]);
if (NIPJNK) contributingNodeIndicesPrCellCorner[3].push_back((*NIPJNK)[5]);
if (NINK ) contributingNodeIndicesPrCellCorner[3].push_back((*NINK )[6]);
// 4 <- NI[5] NINJ[6] NJ[7] PK[0] NIPK[1] NINJPK[2] NJPK[3]
if (IJK ) contributingNodeIndicesPrCellCorner[4].push_back((*IJK )[4]);
if (NI ) contributingNodeIndicesPrCellCorner[4].push_back((*NI )[5]);
if (NINJ ) contributingNodeIndicesPrCellCorner[4].push_back((*NINJ )[6]);
if (NJ ) contributingNodeIndicesPrCellCorner[4].push_back((*NJ )[7]);
if (PK ) contributingNodeIndicesPrCellCorner[4].push_back((*PK )[0]);
if (NIPK ) contributingNodeIndicesPrCellCorner[4].push_back((*NIPK )[1]);
if (NINJPK) contributingNodeIndicesPrCellCorner[4].push_back((*NINJPK)[2]);
if (NJPK ) contributingNodeIndicesPrCellCorner[4].push_back((*NJPK )[3]);
if (IJK ) contributingNodeIndicesPrCellCorner[5].push_back((*IJK )[5]);
if (NJ ) contributingNodeIndicesPrCellCorner[5].push_back((*NJ )[6]);
if (PINJ ) contributingNodeIndicesPrCellCorner[5].push_back((*PINJ )[7]);
if (PI ) contributingNodeIndicesPrCellCorner[5].push_back((*PI )[4]);
if (PK ) contributingNodeIndicesPrCellCorner[5].push_back((*PK )[1]);
if (NJPK ) contributingNodeIndicesPrCellCorner[5].push_back((*NJPK )[2]);
if (PINJPK) contributingNodeIndicesPrCellCorner[5].push_back((*PINJPK)[3]);
if (PIPK ) contributingNodeIndicesPrCellCorner[5].push_back((*PIPK )[0]);
// 6 <- PI[7] PIPJ[4] PJ[5] PK[2] PIPK[3] PIPJPK[0] PJPK[1]
if (IJK ) contributingNodeIndicesPrCellCorner[6].push_back((*IJK )[6]);
if (PI ) contributingNodeIndicesPrCellCorner[6].push_back((*PI )[7]);
if (PIPJ ) contributingNodeIndicesPrCellCorner[6].push_back((*PIPJ )[4]);
if (PJ ) contributingNodeIndicesPrCellCorner[6].push_back((*PJ )[5]);
if (PK ) contributingNodeIndicesPrCellCorner[6].push_back((*PK )[2]);
if (PIPK ) contributingNodeIndicesPrCellCorner[6].push_back((*PIPK )[3]);
if (PIPJPK) contributingNodeIndicesPrCellCorner[6].push_back((*PIPJPK)[0]);
if (PJPK ) contributingNodeIndicesPrCellCorner[6].push_back((*PJPK )[1]);
if (IJK ) contributingNodeIndicesPrCellCorner[7].push_back((*IJK )[7]);
if (PJ ) contributingNodeIndicesPrCellCorner[7].push_back((*PJ )[4]);
if (NIPJ ) contributingNodeIndicesPrCellCorner[7].push_back((*NIPJ )[5]);
if (NI ) contributingNodeIndicesPrCellCorner[7].push_back((*NI )[6]);
if (PK ) contributingNodeIndicesPrCellCorner[7].push_back((*PK )[3]);
if (PJPK ) contributingNodeIndicesPrCellCorner[7].push_back((*PJPK )[0]);
if (NIPJPK) contributingNodeIndicesPrCellCorner[7].push_back((*NIPJPK)[1]);
if (NIPK ) contributingNodeIndicesPrCellCorner[7].push_back((*NIPK )[2]);
// Average the nodes
for (size_t cornIdx = 0; cornIdx < 8; ++cornIdx)
{
estimatedElmCorners[cornIdx] = cvf::Vec3d::ZERO;
size_t contribCount = contributingNodeIndicesPrCellCorner[cornIdx].size();
for (size_t ctnIdx = 0; ctnIdx < contribCount; ++ctnIdx)
{
estimatedElmCorners[cornIdx] += eclNodes[contributingNodeIndicesPrCellCorner[cornIdx][ctnIdx]];
}
estimatedElmCorners[cornIdx] /= contribCount;
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigCaseToCaseCellMapperTools::rotateQuad(cvf::Vec3d quad[4], int idxToNewStart)
{
if (idxToNewStart == 0) return;
cvf::Vec3d tmpQuad[4];
tmpQuad[0] = quad[0];
tmpQuad[1] = quad[1];
tmpQuad[2] = quad[2];
tmpQuad[3] = quad[3];
quad[0] = tmpQuad[idxToNewStart];
++idxToNewStart; if (idxToNewStart > 3) idxToNewStart = 0;
quad[1] = tmpQuad[idxToNewStart];
++idxToNewStart; if (idxToNewStart > 3) idxToNewStart = 0;
quad[2] = tmpQuad[idxToNewStart];
++idxToNewStart; if (idxToNewStart > 3) idxToNewStart = 0;
quad[3] = tmpQuad[idxToNewStart];
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigCaseToCaseCellMapperTools::flipQuadWinding(cvf::Vec3d quad[4])
{
cvf::Vec3d temp = quad[1];
quad[1] = quad[3];
quad[3] = temp;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
int RigCaseToCaseCellMapperTools::quadVxClosestToXYOfPoint( const cvf::Vec3d point, const cvf::Vec3d quad[4])
{
double minSqDist = HUGE_VAL;
int quadVxIdxClosestToPoint = cvf::UNDEFINED_INT;
for (int i = 0; i < 4; ++i)
{
cvf::Vec3d diff = quad[i]- point;
diff[2] = 0.0;
double sqDist = diff.lengthSquared();
if (sqDist < minSqDist)
{
minSqDist = sqDist;
quadVxIdxClosestToPoint = i;
}
}
return quadVxIdxClosestToPoint;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigCaseToCaseCellMapperTools::elementCorners(const RigFemPart* femPart, int elmIdx, cvf::Vec3d elmCorners[8])
{
RigElementType elmType = femPart->elementType(elmIdx);
if (!(elmType == HEX8 || elmType == HEX8P)) return false;
const std::vector<cvf::Vec3f>& nodeCoords = femPart->nodes().coordinates;
const int* cornerIndices = femPart->connectivities(elmIdx);
elmCorners[0] = cvf::Vec3d(nodeCoords[cornerIndices[0]]);
elmCorners[1] = cvf::Vec3d(nodeCoords[cornerIndices[1]]);
elmCorners[2] = cvf::Vec3d(nodeCoords[cornerIndices[2]]);
elmCorners[3] = cvf::Vec3d(nodeCoords[cornerIndices[3]]);
elmCorners[4] = cvf::Vec3d(nodeCoords[cornerIndices[4]]);
elmCorners[5] = cvf::Vec3d(nodeCoords[cornerIndices[5]]);
elmCorners[6] = cvf::Vec3d(nodeCoords[cornerIndices[6]]);
elmCorners[7] = cvf::Vec3d(nodeCoords[cornerIndices[7]]);
return true;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
int RigCaseToCaseCellMapperTools::findMatchingPOSKFaceIdx(const cvf::Vec3d baseCell[8],bool isBaseCellNormalsOutwards, const cvf::Vec3d c2[8])
{
int faceNodeCount;
const int* posKFace = RigFemTypes::localElmNodeIndicesForFace(HEX8, (int)(cvf::StructGridInterface::POS_K), &faceNodeCount);
double sign = isBaseCellNormalsOutwards ? 1.0 : -1.0;
cvf::Vec3d posKnormal = sign*(baseCell[posKFace[2]] - baseCell[posKFace[0]]) ^ (baseCell[posKFace[3]] - baseCell[posKFace[1]]);
posKnormal.normalize();
double minDiff = HUGE_VAL;
int bestFace = -1;
for (int faceIdx = 5; faceIdx >= 0; --faceIdx) // Backwards. might hit earlier more often
{
const int* face = RigFemTypes::localElmNodeIndicesForFace(HEX8, faceIdx, &faceNodeCount);
cvf::Vec3d normal = (c2[face[2]] - c2[face[0]]) ^ (c2[face[3]] - c2[face[1]]);
normal.normalize();
double sqDiff = (posKnormal-normal).lengthSquared();
if (sqDiff < minDiff)
{
minDiff = sqDiff;
bestFace = faceIdx;
if (minDiff < 0.1*0.1) break; // This must be the one. Do not search further
}
}
return bestFace;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigCaseToCaseCellMapperTools::isEclFemCellsMatching(const cvf::Vec3d baseCell[8],
cvf::Vec3d cell[8],
double xyTolerance, double zTolerance)
{
bool isMatching = true;
for (int i = 0; i < 4 ; ++i)
{
cvf::Vec3d diff = cell[i] - baseCell[i];
if (!(fabs(diff.x()) < xyTolerance && fabs(diff.y()) < xyTolerance && fabs(diff.z()) < zTolerance))
{
isMatching = false;
break;
}
}
return isMatching;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigCaseToCaseCellMapperTools::rotateCellTopologicallyToMatchBaseCell(const cvf::Vec3d * baseCell, bool baseCellFaceNormalsIsOutwards, cvf::Vec3d * cell)
{
int femDeepZFaceIdx = findMatchingPOSKFaceIdx(baseCell, baseCellFaceNormalsIsOutwards, cell);
{
cvf::Vec3d tmpFemCorners[8];
tmpFemCorners[0] = cell[0];
tmpFemCorners[1] = cell[1];
tmpFemCorners[2] = cell[2];
tmpFemCorners[3] = cell[3];
tmpFemCorners[4] = cell[4];
tmpFemCorners[5] = cell[5];
tmpFemCorners[6] = cell[6];
tmpFemCorners[7] = cell[7];
int femShallowZFaceIdx = RigFemTypes::oppositeFace(HEX8, femDeepZFaceIdx);
int faceNodeCount;
const int* localElmNodeIndicesForPOSKFace = RigFemTypes::localElmNodeIndicesForFace(HEX8, femDeepZFaceIdx, &faceNodeCount);
const int* localElmNodeIndicesForNEGKFace = RigFemTypes::localElmNodeIndicesForFace(HEX8, femShallowZFaceIdx, &faceNodeCount);
cell[0] = tmpFemCorners[localElmNodeIndicesForNEGKFace[0]];
cell[1] = tmpFemCorners[localElmNodeIndicesForNEGKFace[1]];
cell[2] = tmpFemCorners[localElmNodeIndicesForNEGKFace[2]];
cell[3] = tmpFemCorners[localElmNodeIndicesForNEGKFace[3]];
cell[4] = tmpFemCorners[localElmNodeIndicesForPOSKFace[0]];
cell[5] = tmpFemCorners[localElmNodeIndicesForPOSKFace[1]];
cell[6] = tmpFemCorners[localElmNodeIndicesForPOSKFace[2]];
cell[7] = tmpFemCorners[localElmNodeIndicesForPOSKFace[3]];
}
cvf::Vec3d* femDeepestQuad = &(cell[4]);
cvf::Vec3d* femShallowQuad = &(cell[0]);
// Now the top/bottom have opposite winding. To make the comparisons and index rotations simpler
// flip the winding of the top or bottom face depending on whether the eclipse grid is inside-out
if (baseCellFaceNormalsIsOutwards)
{
flipQuadWinding(femShallowQuad);
}
else
{
flipQuadWinding(femDeepestQuad);
}
// We now need to rotate the fem quads to be alligned with the ecl quads
// Since the start point of the quad always is aligned with the opposite face-quad start
// we can find the rotation for the top, and apply it to both top and bottom
int femQuadStartIdx = quadVxClosestToXYOfPoint(baseCell[0], femShallowQuad);
rotateQuad(femDeepestQuad, femQuadStartIdx);
rotateQuad(femShallowQuad, femQuadStartIdx);
}
#if 0 // Inside Bounding box test
cvf::BoundingBox cellBBox;
for (int i = 0; i < 8 ; ++i) cellBBox.add(cellCorners[i]);
cvf::Vec3d cs = cellBBox.min();
cvf::Vec3d cl = cellBBox.max();
cvf::Vec3d es = elmBBox.min();
cvf::Vec3d el = elmBBox.max();
if ( ( (cs.x() + xyTolerance) >= es.x() && (cl.x() - xyTolerance) <= el.x())
&& ( (cs.y() + xyTolerance) >= es.y() && (cl.y() - xyTolerance) <= el.y())
&& ( (cs.z() + zTolerance ) >= es.z() && (cl.z() - zTolerance ) <= el.z()) )
{
// Cell bb equal or inside Elm bb
isMatching = true;
}
if ( ( (es.x() + xyTolerance) >= cs.x() && (el.x() - xyTolerance) <= cl.x())
&& ( (es.y() + xyTolerance) >= cs.y() && (el.y() - xyTolerance) <= cl.y())
&& ( (es.z() + zTolerance ) >= cs.z() && (el.z() - zTolerance ) <= cl.z()) )
{
// Elm bb equal or inside Cell bb
isMatching = true;
}
#endif
#if 0
{
const std::vector<cvf::Vec3d>& eclNodes = eclGrid->nodes();
const RigCell& cell = eclGrid->cells()[reservoirCellIndex];
const caf::SizeTArray8& cornerIndices = cell.cornerIndices();
int faceNodeCount;
const int* localElmNodeIndicesForTopZFace = RigFemTypes::localElmNodeIndicesForFace(HEX8, 4, &faceNodeCount);
const int* localElmNodeIndicesForBotZFace = RigFemTypes::localElmNodeIndicesForFace(HEX8, 5, &faceNodeCount);
eclDeepestQuad[0] = eclNodes[cornerIndices[localElmNodeIndicesForTopZFace[0]]];
eclDeepestQuad[1] = eclNodes[cornerIndices[localElmNodeIndicesForTopZFace[1]]];
eclDeepestQuad[2] = eclNodes[cornerIndices[localElmNodeIndicesForTopZFace[2]]];
eclDeepestQuad[3] = eclNodes[cornerIndices[localElmNodeIndicesForTopZFace[3]]];
eclShallowQuad[0] = eclNodes[cornerIndices[localElmNodeIndicesForBotZFace[0]]];
eclShallowQuad[1] = eclNodes[cornerIndices[localElmNodeIndicesForBotZFace[1]]];
eclShallowQuad[2] = eclNodes[cornerIndices[localElmNodeIndicesForBotZFace[2]]];
eclShallowQuad[3] = eclNodes[cornerIndices[localElmNodeIndicesForBotZFace[3]]];
}
#endif
#if 0
// First search K=1 diagonally for a seed cell; A cell without collapsings, and without faults
size_t minIJCount = masterEclGrid->cellCountI();
if (minIJCount > masterEclGrid->cellCountJ())
minIJCount = masterEclGrid->cellCountJ();
for (size_t ij = 0; ij < minIJCount; ++ij )
{
size_t localCellIdx = masterEclGrid->cellIndexFromIJK(ij, ij, 0);
size_t reservoirCellIdx = masterEclGrid->reservoirCellIndex(localCellIdx);
cvf::Vec3d vertices[8];
masterEclGrid->cellCornerVertices(localCellIdx, vertices);
if (!isCellNormal(vertices))
continue;
const RigFault* fault = masterEclGrid->findFaultFromCellIndexAndCellFace(reservoirCellIdx, cvf::StructGridInterface::POS_I);
}
#endif
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigCaseToCaseCellMapperTools::calculateCellCenter(cvf::Vec3d elmCorners[8])
{
cvf::Vec3d avg(cvf::Vec3d::ZERO);
size_t i;
for (i = 0; i < 8; i++)
{
avg += elmCorners[i];
}
avg /= 8.0;
return avg;
}

54
ApplicationCode/ReservoirDataModel/RigCaseToCaseCellMapperTools.h Normal file
View File

@ -0,0 +1,54 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2015- Statoil ASA
// Copyright (C) 2015- Ceetron Solutions AS
//
// 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 "cvfObject.h"
#include "cvfMath.h"
#include "cvfVector3.h"
#include <vector>
class RigMainGrid;
class RigFemPart;
//==================================================================================================
///
//==================================================================================================
class RigCaseToCaseCellMapperTools
{
public:
static void estimatedFemCellFromEclCell(const RigMainGrid* eclGrid, size_t reservoirCellIndex, cvf::Vec3d estimatedElmCorners[8]);
static void rotateCellTopologicallyToMatchBaseCell(const cvf::Vec3d * baseCell, bool baseCellFaceNormalsIsOutwards, cvf::Vec3d * cell);
static cvf::Vec3d calculateCellCenter(cvf::Vec3d elmCorners[8]);
static bool elementCorners(const RigFemPart* femPart, int elmIdx, cvf::Vec3d elmCorners[8]);
static bool isEclFemCellsMatching(const cvf::Vec3d baseCell[8], cvf::Vec3d cell[8], double xyTolerance, double zTolerance);
private:
static void rotateQuad(cvf::Vec3d quad[4], int idxToNewStart);
static void flipQuadWinding(cvf::Vec3d quad[4]);
static int quadVxClosestToXYOfPoint(const cvf::Vec3d point, const cvf::Vec3d quad[4]);
static int findMatchingPOSKFaceIdx(const cvf::Vec3d baseCell[8], bool isBaseCellNormalsOutwards, const cvf::Vec3d c2[8]);
};

295
ApplicationCode/ReservoirDataModel/RigCaseToCaseRangeFilterMapper.cpp Normal file
View File

@ -0,0 +1,295 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2015- Statoil ASA
// Copyright (C) 2015- Ceetron Solutions AS
//
// 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 "RigCaseToCaseRangeFilterMapper.h"
#include "RigCaseToCaseCellMapper.h"
#include "RigCaseToCaseCellMapperTools.h"
#include "RigFemPart.h"
#include "RigMainGrid.h"
#include "RigFemPartGrid.h"
#include "RimCellRangeFilter.h"
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigCaseToCaseRangeFilterMapper::convertRangeFilterEclToFem(RimCellRangeFilter* srcFilter, const RigMainGrid* srcEclGrid, RimCellRangeFilter* dstFilter, const RigFemPart* dstFemPart)
{
CVF_ASSERT(srcFilter && srcEclGrid && dstFilter && dstFemPart);
CVF_ASSERT(srcFilter->gridIndex() == 0); // LGR not supported yet
struct RangeFilterCorner { RangeFilterCorner() : isExactMatch(false){} cvf::Vec3st ijk; bool isExactMatch; };
RangeFilterCorner rangeFilterMatches[8];
size_t srcStartI = srcFilter->startIndexI() - 1;
size_t srcStartJ = srcFilter->startIndexJ() - 1;
size_t srcStartK = srcFilter->startIndexK() - 1;
size_t srcEndI = srcStartI + srcFilter->cellCountI();
size_t srcEndJ = srcStartJ + srcFilter->cellCountJ();
size_t srcEndK = srcStartK + srcFilter->cellCountK();
cvf::Vec3st srcRangeCube[8];
srcRangeCube[0] = cvf::Vec3st(srcStartI, srcStartJ, srcStartK);
srcRangeCube[1] = cvf::Vec3st(srcEndI, srcStartJ, srcStartK);
srcRangeCube[2] = cvf::Vec3st(srcEndI, srcEndJ, srcStartK);
srcRangeCube[3] = cvf::Vec3st(srcStartI, srcEndJ, srcStartK);
srcRangeCube[4] = cvf::Vec3st(srcStartI, srcStartJ, srcEndK);
srcRangeCube[5] = cvf::Vec3st(srcEndI, srcStartJ, srcEndK);
srcRangeCube[6] = cvf::Vec3st(srcEndI, srcEndJ, srcEndK);
srcRangeCube[7] = cvf::Vec3st(srcStartI, srcEndJ, srcEndK);
size_t dstStartI = cvf::UNDEFINED_SIZE_T;
size_t dstStartJ = cvf::UNDEFINED_SIZE_T;
size_t dstStartK = cvf::UNDEFINED_SIZE_T;
size_t dstEndI = cvf::UNDEFINED_SIZE_T;
size_t dstEndJ = cvf::UNDEFINED_SIZE_T;
size_t dstEndK = cvf::UNDEFINED_SIZE_T;
bool foundExactMatch = false;
int cornerIdx = 0;
int diagIdx = 6;// Index to diagonal corner
for (cornerIdx = 0; cornerIdx < 4; ++cornerIdx)
{
diagIdx = (cornerIdx < 2) ? cornerIdx + 6 : cornerIdx + 2;
rangeFilterMatches[cornerIdx].isExactMatch = findBestFemCellFromEclCell(srcEclGrid,
srcRangeCube[cornerIdx][0],
srcRangeCube[cornerIdx][1],
srcRangeCube[cornerIdx][2],
dstFemPart,
&(rangeFilterMatches[cornerIdx].ijk[0]),
&(rangeFilterMatches[cornerIdx].ijk[1]),
&(rangeFilterMatches[cornerIdx].ijk[2]));
rangeFilterMatches[diagIdx].isExactMatch = findBestFemCellFromEclCell(srcEclGrid,
srcRangeCube[diagIdx][0],
srcRangeCube[diagIdx][1],
srcRangeCube[diagIdx][2],
dstFemPart,
&(rangeFilterMatches[diagIdx].ijk[0]),
&(rangeFilterMatches[diagIdx].ijk[1]),
&(rangeFilterMatches[diagIdx].ijk[2]));
if (rangeFilterMatches[cornerIdx].isExactMatch && rangeFilterMatches[diagIdx].isExactMatch)
{
foundExactMatch = true;
break;
}
}
// Get the start and end IJK from the matched corners
if (foundExactMatch)
{
// Populate dst range filter from the diagonal that matches exact
dstStartI = CVF_MIN(rangeFilterMatches[cornerIdx].ijk[0], rangeFilterMatches[diagIdx].ijk[0]);
dstStartJ = CVF_MIN(rangeFilterMatches[cornerIdx].ijk[1], rangeFilterMatches[diagIdx].ijk[1]);
dstStartK = CVF_MIN(rangeFilterMatches[cornerIdx].ijk[2], rangeFilterMatches[diagIdx].ijk[2]);
dstEndI = CVF_MAX(rangeFilterMatches[cornerIdx].ijk[0], rangeFilterMatches[diagIdx].ijk[0]);
dstEndJ = CVF_MAX(rangeFilterMatches[cornerIdx].ijk[1], rangeFilterMatches[diagIdx].ijk[1]);
dstEndK = CVF_MAX(rangeFilterMatches[cornerIdx].ijk[2], rangeFilterMatches[diagIdx].ijk[2]);
}
else
{
// Todo: be even smarter, and use possible matching corners to add up an as best solution as possible.
// For now we just take the first diagonal.
dstStartI = rangeFilterMatches[0].ijk[0];
dstStartJ = rangeFilterMatches[0].ijk[1];
dstStartK = rangeFilterMatches[0].ijk[2];
dstEndI = rangeFilterMatches[6].ijk[0];
dstEndJ = rangeFilterMatches[6].ijk[1];
dstEndK = rangeFilterMatches[6].ijk[2];
}
// Populate the dst range filter with new data
if ((dstStartI != cvf::UNDEFINED_SIZE_T && dstStartJ != cvf::UNDEFINED_SIZE_T && dstStartK != cvf::UNDEFINED_SIZE_T)
&& (dstEndI != cvf::UNDEFINED_SIZE_T && dstEndJ != cvf::UNDEFINED_SIZE_T && dstEndK != cvf::UNDEFINED_SIZE_T))
{
dstFilter->startIndexJ = static_cast<int>(dstStartI + 1);
dstFilter->startIndexK = static_cast<int>(dstStartJ + 1);
dstFilter->startIndexI = static_cast<int>(dstStartK + 1);
dstFilter->cellCountI = static_cast<int>(dstEndI - (dstStartI-1));
dstFilter->cellCountJ = static_cast<int>(dstEndJ - (dstStartJ-1));
dstFilter->cellCountK = static_cast<int>(dstEndK - (dstStartK-1));
}
else
{
dstFilter->startIndexI = 1;
dstFilter->startIndexJ = 1;
dstFilter->startIndexK = 1;
dstFilter->cellCountI = 0;
dstFilter->cellCountJ = 0;
dstFilter->cellCountK = 0;
dstFilter->computeAndSetValidValues();
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigCaseToCaseRangeFilterMapper::findBestFemCellFromEclCell(const RigMainGrid* masterEclGrid, size_t ei, size_t ej, size_t ek, const RigFemPart* dependentFemPart, size_t* fi, size_t * fj, size_t* fk)
{
// Find tolerance
double cellSizeI, cellSizeJ, cellSizeK;
masterEclGrid->characteristicCellSizes(&cellSizeI, &cellSizeJ, &cellSizeK);
double xyTolerance = cellSizeI* 0.4;
double zTolerance = cellSizeK* 0.4;
bool isEclFaceNormalsOutwards = masterEclGrid->isFaceNormalsOutwards();
size_t cellIdx = masterEclGrid->cellIndexFromIJK(ei, ej, ek);
cvf::Vec3d geoMechConvertedEclCell[8];
RigCaseToCaseCellMapperTools::estimatedFemCellFromEclCell(masterEclGrid, cellIdx, geoMechConvertedEclCell);
cvf::BoundingBox elmBBox;
for (int i = 0; i < 8 ; ++i) elmBBox.add(geoMechConvertedEclCell[i]);
std::vector<size_t> closeElements;
dependentFemPart->findIntersectingCells(elmBBox, &closeElements);
cvf::Vec3d elmCorners[8];
int elmIdxToBestMatch = -1;
double sqDistToClosestElmCenter = HUGE_VAL;
cvf::Vec3d convEclCellCenter = RigCaseToCaseCellMapperTools::calculateCellCenter(geoMechConvertedEclCell);
bool foundExactMatch = false;
for (size_t ccIdx = 0; ccIdx < closeElements.size(); ++ccIdx)
{
int elmIdx = static_cast<int>(closeElements[ccIdx]);
RigCaseToCaseCellMapperTools::elementCorners(dependentFemPart, elmIdx, elmCorners);
cvf::Vec3d cellCenter = RigCaseToCaseCellMapperTools::calculateCellCenter(elmCorners);
double sqDist = (cellCenter - convEclCellCenter).lengthSquared();
if (sqDist < sqDistToClosestElmCenter)
{
elmIdxToBestMatch = elmIdx;
sqDistToClosestElmCenter = sqDist;
}
RigCaseToCaseCellMapperTools::rotateCellTopologicallyToMatchBaseCell(geoMechConvertedEclCell, isEclFaceNormalsOutwards, elmCorners);
foundExactMatch = RigCaseToCaseCellMapperTools::isEclFemCellsMatching(geoMechConvertedEclCell, elmCorners,
xyTolerance, zTolerance);
if (foundExactMatch)
{
elmIdxToBestMatch = elmIdx;
break;
}
}
if (elmIdxToBestMatch != -1)
{
dependentFemPart->structGrid()->ijkFromCellIndex(elmIdxToBestMatch, fi, fj, fk);
}
else
{
(*fi) = cvf::UNDEFINED_SIZE_T;
(*fj) = cvf::UNDEFINED_SIZE_T;
(*fk) = cvf::UNDEFINED_SIZE_T;
}
return foundExactMatch;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigCaseToCaseRangeFilterMapper::findBestEclCellFromFemCell(const RigFemPart* dependentFemPart, size_t fi, size_t fj, size_t fk, const RigMainGrid* masterEclGrid, size_t* ei, size_t* ej, size_t* ek)
{
// Find tolerance
double cellSizeI, cellSizeJ, cellSizeK;
masterEclGrid->characteristicCellSizes(&cellSizeI, &cellSizeJ, &cellSizeK);
double xyTolerance = cellSizeI* 0.4;
double zTolerance = cellSizeK* 0.4;
bool isEclFaceNormalsOutwards = masterEclGrid->isFaceNormalsOutwards();
int elementIdx = static_cast<int>(dependentFemPart->structGrid()->cellIndexFromIJK(fi, fj, fk));
cvf::Vec3d elmCorners[8];
RigCaseToCaseCellMapperTools::elementCorners(dependentFemPart, elementIdx, elmCorners);
cvf::BoundingBox elmBBox;
for (int i = 0; i < 8 ; ++i) elmBBox.add(elmCorners[i]);
std::vector<size_t> closeCells;
masterEclGrid->findIntersectingCells(elmBBox, &closeCells); // This might actually miss the exact one, but we have no other alternative yet.
size_t globCellIdxToBestMatch = cvf::UNDEFINED_SIZE_T;
double sqDistToClosestCellCenter = HUGE_VAL;
cvf::Vec3d elmCenter = RigCaseToCaseCellMapperTools::calculateCellCenter(elmCorners);
bool foundExactMatch = false;
for (size_t ccIdx = 0; ccIdx < closeCells.size(); ++ccIdx)
{
size_t cellIdx = closeCells[ccIdx];
cvf::Vec3d geoMechConvertedEclCell[8];
RigCaseToCaseCellMapperTools::estimatedFemCellFromEclCell(masterEclGrid, cellIdx, geoMechConvertedEclCell);
cvf::Vec3d cellCenter = RigCaseToCaseCellMapperTools::calculateCellCenter(geoMechConvertedEclCell);
double sqDist = (cellCenter - elmCenter).lengthSquared();
if (sqDist < sqDistToClosestCellCenter)
{
globCellIdxToBestMatch = cellIdx;
sqDistToClosestCellCenter = sqDist;
}
RigCaseToCaseCellMapperTools::rotateCellTopologicallyToMatchBaseCell(geoMechConvertedEclCell, isEclFaceNormalsOutwards, elmCorners);
foundExactMatch = RigCaseToCaseCellMapperTools::isEclFemCellsMatching(geoMechConvertedEclCell, elmCorners,
xyTolerance, zTolerance);
if (foundExactMatch)
{
globCellIdxToBestMatch = cellIdx;
break;
}
}
if (globCellIdxToBestMatch != cvf::UNDEFINED_SIZE_T)
{
masterEclGrid->ijkFromCellIndex(globCellIdxToBestMatch, ei, ej, ek);
}
else
{
(*ei) = cvf::UNDEFINED_SIZE_T;
(*ej) = cvf::UNDEFINED_SIZE_T;
(*ek) = cvf::UNDEFINED_SIZE_T;
}
return foundExactMatch;
}

39
ApplicationCode/ReservoirDataModel/RigCaseToCaseRangeFilterMapper.h Normal file
View File

@ -0,0 +1,39 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2015- Statoil ASA
// Copyright (C) 2015- Ceetron Solutions AS
//
// 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
class RimCellRangeFilter;
class RigMainGrid;
class RigFemPart;
class RigCaseToCaseRangeFilterMapper
{
public:
static void convertRangeFilterEclToFem(RimCellRangeFilter* srcFilter, const RigMainGrid* srcEclGrid,
RimCellRangeFilter* dstFilter, const RigFemPart* dstFemPart);
private:
static bool findBestFemCellFromEclCell(const RigMainGrid* masterEclGrid, size_t ei, size_t ej, size_t ek,
const RigFemPart* dependentFemPart, size_t* fi, size_t * fj, size_t* fk);
static bool findBestEclCellFromFemCell(const RigFemPart* dependentFemPart, size_t fi, size_t fj, size_t fk,
const RigMainGrid* masterEclGrid, size_t* ei, size_t* ej, size_t* ek);
};