Export Completion: Fix missing segments and connections for LGR grids

* Compute characteristic cell size based on active cells
* Compute cell face normal based on a suitable active cells
Using all cells as basis for face normal is fragile. Some models have highly distorted/twisted cells, but all active cells should be geometrically more stable.
This commit is contained in:
Magne Sjaastad
2023-02-15 15:25:38 +01:00
committed by GitHub
parent a7acbe2c86
commit 20439e1da9
5 changed files with 163 additions and 98 deletions

View File

@@ -38,6 +38,8 @@
#include "cvfBase.h"
#include "cvfBoundingBox.h"
#include <algorithm>
namespace caf
{
template <>
@@ -319,99 +321,13 @@ void StructGridInterface::characteristicCellSizes( double* iSize, double* jSize,
{
CVF_ASSERT( iSize && jSize && kSize );
if ( m_characteristicCellSizeI == cvf::UNDEFINED_DOUBLE || m_characteristicCellSizeJ == cvf::UNDEFINED_DOUBLE ||
m_characteristicCellSizeK == cvf::UNDEFINED_DOUBLE )
if ( !hasValidCharacteristicCellSizes() )
{
ubyte faceConnPosI[4];
cellFaceVertexIndices( StructGridInterface::POS_I, faceConnPosI );
std::vector<size_t> reservoirCellIndices;
reservoirCellIndices.resize( cellCountI() * cellCountJ() * cellCountK() );
std::iota( reservoirCellIndices.begin(), reservoirCellIndices.end(), 0 );
ubyte faceConnNegI[4];
cellFaceVertexIndices( StructGridInterface::NEG_I, faceConnNegI );
ubyte faceConnPosJ[4];
cellFaceVertexIndices( StructGridInterface::POS_J, faceConnPosJ );
ubyte faceConnNegJ[4];
cellFaceVertexIndices( StructGridInterface::NEG_J, faceConnNegJ );
ubyte faceConnPosK[4];
cellFaceVertexIndices( StructGridInterface::POS_K, faceConnPosK );
ubyte faceConnNegK[4];
cellFaceVertexIndices( StructGridInterface::NEG_K, faceConnNegK );
double iLengthAccumulated = 0.0;
double jLengthAccumulated = 0.0;
double kLengthAccumulated = 0.0;
cvf::Vec3d cornerVerts[8];
size_t cellCount = 0;
size_t k;
for ( k = 0; k < cellCountK(); k++ )
{
size_t j;
for ( j = 0; j < cellCountJ(); j++ )
{
size_t i;
for ( i = 0; i < cellCountI(); i += 10 ) // NB! Evaluate every n-th cell
{
if ( isCellValid( i, j, k ) )
{
size_t cellIndex = cellIndexFromIJK( i, j, k );
cellCornerVertices( cellIndex, cornerVerts );
cvf::BoundingBox bb;
for ( const auto& v : cornerVerts )
{
bb.add( v );
}
// Exclude cells with very small volumes
const double tolerance = 0.2;
if ( bb.extent().z() < tolerance ) continue;
iLengthAccumulated +=
( cornerVerts[faceConnPosI[0]] - cornerVerts[faceConnNegI[0]] ).lengthSquared();
iLengthAccumulated +=
( cornerVerts[faceConnPosI[1]] - cornerVerts[faceConnNegI[3]] ).lengthSquared();
iLengthAccumulated +=
( cornerVerts[faceConnPosI[2]] - cornerVerts[faceConnNegI[2]] ).lengthSquared();
iLengthAccumulated +=
( cornerVerts[faceConnPosI[3]] - cornerVerts[faceConnNegI[1]] ).lengthSquared();
jLengthAccumulated +=
( cornerVerts[faceConnPosJ[0]] - cornerVerts[faceConnNegJ[0]] ).lengthSquared();
jLengthAccumulated +=
( cornerVerts[faceConnPosJ[1]] - cornerVerts[faceConnNegJ[3]] ).lengthSquared();
jLengthAccumulated +=
( cornerVerts[faceConnPosJ[2]] - cornerVerts[faceConnNegJ[2]] ).lengthSquared();
jLengthAccumulated +=
( cornerVerts[faceConnPosJ[3]] - cornerVerts[faceConnNegJ[1]] ).lengthSquared();
kLengthAccumulated +=
( cornerVerts[faceConnPosK[0]] - cornerVerts[faceConnNegK[0]] ).lengthSquared();
kLengthAccumulated +=
( cornerVerts[faceConnPosK[1]] - cornerVerts[faceConnNegK[3]] ).lengthSquared();
kLengthAccumulated +=
( cornerVerts[faceConnPosK[2]] - cornerVerts[faceConnNegK[2]] ).lengthSquared();
kLengthAccumulated +=
( cornerVerts[faceConnPosK[3]] - cornerVerts[faceConnNegK[1]] ).lengthSquared();
cellCount++;
}
}
}
}
double divisor = cellCount * 4.0;
if ( divisor > 0.0 )
{
m_characteristicCellSizeI = cvf::Math::sqrt( iLengthAccumulated / divisor );
m_characteristicCellSizeJ = cvf::Math::sqrt( jLengthAccumulated / divisor );
m_characteristicCellSizeK = cvf::Math::sqrt( kLengthAccumulated / divisor );
}
computeCharacteristicCellSize( reservoirCellIndices );
}
*iSize = m_characteristicCellSizeI;
@@ -419,4 +335,100 @@ void StructGridInterface::characteristicCellSizes( double* iSize, double* jSize,
*kSize = m_characteristicCellSizeK;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool StructGridInterface::hasValidCharacteristicCellSizes() const
{
if ( m_characteristicCellSizeI == cvf::UNDEFINED_DOUBLE || m_characteristicCellSizeJ == cvf::UNDEFINED_DOUBLE ||
m_characteristicCellSizeK == cvf::UNDEFINED_DOUBLE )
return false;
return true;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void StructGridInterface::computeCharacteristicCellSize( const std::vector<size_t>& globalCellIndices ) const
{
ubyte faceConnPosI[4];
cellFaceVertexIndices( StructGridInterface::POS_I, faceConnPosI );
ubyte faceConnNegI[4];
cellFaceVertexIndices( StructGridInterface::NEG_I, faceConnNegI );
ubyte faceConnPosJ[4];
cellFaceVertexIndices( StructGridInterface::POS_J, faceConnPosJ );
ubyte faceConnNegJ[4];
cellFaceVertexIndices( StructGridInterface::NEG_J, faceConnNegJ );
ubyte faceConnPosK[4];
cellFaceVertexIndices( StructGridInterface::POS_K, faceConnPosK );
ubyte faceConnNegK[4];
cellFaceVertexIndices( StructGridInterface::NEG_K, faceConnNegK );
double iLengthAccumulated = 0.0;
double jLengthAccumulated = 0.0;
double kLengthAccumulated = 0.0;
cvf::Vec3d cornerVerts[8];
size_t evaluatedCellCount = 0;
// Evaluate N-th cells, compute the stride between each index
size_t stride = std::max( size_t( 1 ), globalCellIndices.size() / 100 );
size_t i, j, k = 0;
size_t index = 0;
while ( index < globalCellIndices.size() - 1 )
{
size_t cellIndex = globalCellIndices[index];
ijkFromCellIndex( cellIndex, &i, &j, &k );
if ( isCellValid( i, j, k ) )
{
cellCornerVertices( cellIndex, cornerVerts );
cvf::BoundingBox bb;
for ( const auto& v : cornerVerts )
{
bb.add( v );
}
// Exclude cells with very small volumes
const double tolerance = 0.2;
if ( bb.extent().z() < tolerance ) continue;
iLengthAccumulated += ( cornerVerts[faceConnPosI[0]] - cornerVerts[faceConnNegI[0]] ).lengthSquared();
iLengthAccumulated += ( cornerVerts[faceConnPosI[1]] - cornerVerts[faceConnNegI[3]] ).lengthSquared();
iLengthAccumulated += ( cornerVerts[faceConnPosI[2]] - cornerVerts[faceConnNegI[2]] ).lengthSquared();
iLengthAccumulated += ( cornerVerts[faceConnPosI[3]] - cornerVerts[faceConnNegI[1]] ).lengthSquared();
jLengthAccumulated += ( cornerVerts[faceConnPosJ[0]] - cornerVerts[faceConnNegJ[0]] ).lengthSquared();
jLengthAccumulated += ( cornerVerts[faceConnPosJ[1]] - cornerVerts[faceConnNegJ[3]] ).lengthSquared();
jLengthAccumulated += ( cornerVerts[faceConnPosJ[2]] - cornerVerts[faceConnNegJ[2]] ).lengthSquared();
jLengthAccumulated += ( cornerVerts[faceConnPosJ[3]] - cornerVerts[faceConnNegJ[1]] ).lengthSquared();
kLengthAccumulated += ( cornerVerts[faceConnPosK[0]] - cornerVerts[faceConnNegK[0]] ).lengthSquared();
kLengthAccumulated += ( cornerVerts[faceConnPosK[1]] - cornerVerts[faceConnNegK[3]] ).lengthSquared();
kLengthAccumulated += ( cornerVerts[faceConnPosK[2]] - cornerVerts[faceConnNegK[2]] ).lengthSquared();
kLengthAccumulated += ( cornerVerts[faceConnPosK[3]] - cornerVerts[faceConnNegK[1]] ).lengthSquared();
evaluatedCellCount++;
}
index += stride;
}
double divisor = evaluatedCellCount * 4.0;
if ( divisor > 0.0 )
{
m_characteristicCellSizeI = cvf::Math::sqrt( iLengthAccumulated / divisor );
m_characteristicCellSizeJ = cvf::Math::sqrt( jLengthAccumulated / divisor );
m_characteristicCellSizeK = cvf::Math::sqrt( kLengthAccumulated / divisor );
}
}
} // namespace cvf

View File

@@ -91,6 +91,9 @@ public:
virtual cvf::Vec3d maxCoordinate() const = 0;
void characteristicCellSizes( double* iSize, double* jSize, double* kSize ) const;
bool hasValidCharacteristicCellSizes() const;
void computeCharacteristicCellSize( const std::vector<size_t>& globalCellIndices ) const;
virtual cvf::Vec3d displayModelOffset() const;
virtual bool cellIJKNeighbor( size_t i, size_t j, size_t k, FaceType face, size_t* neighborCellIndex ) const = 0;