mirror of
https://github.com/OPM/ResInsight.git
synced 2025-01-01 03:37:15 -06:00
481 lines
19 KiB
C++
481 lines
19 KiB
C++
/////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Copyright (C) 2011- Statoil ASA
|
|
// Copyright (C) 2013- Ceetron Solutions AS
|
|
// Copyright (C) 2011-2012 Ceetron 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 "RigCell.h"
|
|
#include "RigCellGeometryTools.h"
|
|
#include "RigMainGrid.h"
|
|
#include "cvfPlane.h"
|
|
#include "cvfRay.h"
|
|
|
|
#include <cmath>
|
|
|
|
static size_t undefinedCornersArray[8] = { cvf::UNDEFINED_SIZE_T,
|
|
cvf::UNDEFINED_SIZE_T,
|
|
cvf::UNDEFINED_SIZE_T,
|
|
cvf::UNDEFINED_SIZE_T,
|
|
cvf::UNDEFINED_SIZE_T,
|
|
cvf::UNDEFINED_SIZE_T,
|
|
cvf::UNDEFINED_SIZE_T,
|
|
cvf::UNDEFINED_SIZE_T };
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
RigCell::RigCell()
|
|
: m_gridLocalCellIndex( cvf::UNDEFINED_SIZE_T )
|
|
, m_hostGrid( nullptr )
|
|
, m_subGrid( nullptr )
|
|
, m_parentCellIndex( cvf::UNDEFINED_SIZE_T )
|
|
, m_mainGridCellIndex( cvf::UNDEFINED_SIZE_T )
|
|
, m_coarseningBoxIndex( cvf::UNDEFINED_SIZE_T )
|
|
, m_isInvalid( false )
|
|
{
|
|
memcpy( m_cornerIndices.data(), undefinedCornersArray, 8 * sizeof( size_t ) );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
RigCell::~RigCell()
|
|
{
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
cvf::Vec3d RigCell::center() const
|
|
{
|
|
cvf::Vec3d avg( cvf::Vec3d::ZERO );
|
|
|
|
size_t i;
|
|
for ( i = 0; i < 8; i++ )
|
|
{
|
|
avg += m_hostGrid->mainGrid()->nodes()[m_cornerIndices[i]];
|
|
}
|
|
|
|
avg /= 8.0;
|
|
|
|
return avg;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// Get the coordinates of the 4 corners of the given face
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::array<cvf::Vec3d, 4> RigCell::faceCorners( cvf::StructGridInterface::FaceType face ) const
|
|
{
|
|
std::array<cvf::Vec3d, 4> corners;
|
|
cvf::ubyte faceVertexIndices[4];
|
|
cvf::StructGridInterface::cellFaceVertexIndices( face, faceVertexIndices );
|
|
|
|
const std::vector<cvf::Vec3d>& nodeCoords = m_hostGrid->mainGrid()->nodes();
|
|
|
|
for ( size_t i = 0; i < 4; i++ )
|
|
{
|
|
corners[i] = nodeCoords[m_cornerIndices[faceVertexIndices[i]]];
|
|
}
|
|
|
|
return corners;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool isNear( const cvf::Vec3d& p1, const cvf::Vec3d& p2, double tolerance )
|
|
{
|
|
return cvf::Math::abs( p1[0] - p2[0] ) < tolerance && cvf::Math::abs( p1[1] - p2[1] ) < tolerance &&
|
|
cvf::Math::abs( p1[2] - p2[2] ) < tolerance;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool RigCell::isLongPyramidCell( double maxHeightFactor, double nodeNearTolerance ) const
|
|
{
|
|
cvf::ubyte faceVertexIndices[4];
|
|
double squaredMaxHeightFactor = maxHeightFactor * maxHeightFactor;
|
|
|
|
const std::vector<cvf::Vec3d>& nodes = m_hostGrid->mainGrid()->nodes();
|
|
|
|
int face;
|
|
for ( face = 0; face < 6; ++face )
|
|
{
|
|
cvf::StructGridInterface::cellFaceVertexIndices( static_cast<cvf::StructGridInterface::FaceType>( face ), faceVertexIndices );
|
|
int zeroLengthEdgeCount = 0;
|
|
|
|
const cvf::Vec3d& c0 = nodes[m_cornerIndices[faceVertexIndices[0]]];
|
|
const cvf::Vec3d& c1 = nodes[m_cornerIndices[faceVertexIndices[1]]];
|
|
const cvf::Vec3d& c2 = nodes[m_cornerIndices[faceVertexIndices[2]]];
|
|
const cvf::Vec3d& c3 = nodes[m_cornerIndices[faceVertexIndices[3]]];
|
|
|
|
if ( isNear( c0, c1, nodeNearTolerance ) )
|
|
{
|
|
++zeroLengthEdgeCount;
|
|
}
|
|
if ( isNear( c1, c2, nodeNearTolerance ) )
|
|
{
|
|
++zeroLengthEdgeCount;
|
|
}
|
|
if ( isNear( c2, c3, nodeNearTolerance ) )
|
|
{
|
|
++zeroLengthEdgeCount;
|
|
}
|
|
|
|
if ( zeroLengthEdgeCount == 3 )
|
|
{
|
|
return true;
|
|
|
|
#if 0 // More advanced checks turned off since the start. Why did I do that ?
|
|
// Collapse of a complete face is detected. This is possibly the top of a pyramid
|
|
|
|
// "face" has the index to the collapsed face. We need the size of the opposite face
|
|
// to compare it with the pyramid "roof" length.
|
|
|
|
cvf::StructGridInterface::FaceType oppositeFace = cvf::StructGridInterface::POS_I;
|
|
switch (face)
|
|
{
|
|
case cvf::StructGridInterface::POS_I:
|
|
oppositeFace = cvf::StructGridInterface::NEG_I;
|
|
break;
|
|
case cvf::StructGridInterface::POS_J:
|
|
oppositeFace = cvf::StructGridInterface::NEG_J;
|
|
break;
|
|
case cvf::StructGridInterface::POS_K:
|
|
oppositeFace = cvf::StructGridInterface::NEG_K;
|
|
break;
|
|
case cvf::StructGridInterface::NEG_I:
|
|
oppositeFace = cvf::StructGridInterface::POS_I;
|
|
break;
|
|
case cvf::StructGridInterface::NEG_J:
|
|
oppositeFace = cvf::StructGridInterface::POS_J;
|
|
break;
|
|
case cvf::StructGridInterface::NEG_K:
|
|
oppositeFace = cvf::StructGridInterface::POS_K;
|
|
break;
|
|
default:
|
|
CVF_ASSERT(false);
|
|
break;
|
|
}
|
|
|
|
cvf::StructGridInterface::cellFaceVertexIndices(oppositeFace, faceVertexIndices);
|
|
|
|
|
|
const cvf::Vec3d& c0opp = nodes[m_cornerIndices[faceVertexIndices[0]]];
|
|
const cvf::Vec3d& c1opp = nodes[m_cornerIndices[faceVertexIndices[1]]];
|
|
const cvf::Vec3d& c2opp = nodes[m_cornerIndices[faceVertexIndices[2]]];
|
|
const cvf::Vec3d& c3opp = nodes[m_cornerIndices[faceVertexIndices[3]]];
|
|
|
|
// Check if any of the opposite face vertexes are also degenerated to the pyramid top
|
|
|
|
int okVertexCount = 0;
|
|
cvf::Vec3d okVxs[4];
|
|
if (!isNear(c0opp, c0, nodeNearTolerance)) { okVxs[okVertexCount] = c0opp; ++okVertexCount; }
|
|
if (!isNear(c1opp, c0, nodeNearTolerance)) { okVxs[okVertexCount] = c1opp; ++okVertexCount; }
|
|
if (!isNear(c2opp, c0, nodeNearTolerance)) { okVxs[okVertexCount] = c2opp; ++okVertexCount; }
|
|
if (!isNear(c3opp, c0, nodeNearTolerance)) { okVxs[okVertexCount] = c3opp; ++okVertexCount; }
|
|
|
|
if (okVertexCount < 2)
|
|
{
|
|
return true;
|
|
}
|
|
else
|
|
{
|
|
// Use the good vertices to calculate a face size that can be compared to the pyramid height:
|
|
double typicalSquaredEdgeLength = 0;
|
|
for (int i = 1; i < okVertexCount; ++i)
|
|
{
|
|
typicalSquaredEdgeLength += (okVxs[i-1] - okVxs[i]).lengthSquared();
|
|
}
|
|
typicalSquaredEdgeLength /= okVertexCount;
|
|
double pyramidHeightSquared = (okVxs[0] - c0).lengthSquared();
|
|
|
|
if (pyramidHeightSquared > squaredMaxHeightFactor*typicalSquaredEdgeLength)
|
|
{
|
|
return true;
|
|
}
|
|
}
|
|
#endif
|
|
}
|
|
|
|
// Check the ratio of the length of opposite edges.
|
|
// both ratios have to be above threshold to detect a pyramid-ish cell
|
|
// Only test this if we have all nonzero edge lengths.
|
|
else if ( zeroLengthEdgeCount == 0 ) // If the four first faces are ok, the two last must be as well
|
|
{
|
|
double e0SquareLength = ( c1 - c0 ).lengthSquared();
|
|
double e2SquareLength = ( c3 - c2 ).lengthSquared();
|
|
if ( e0SquareLength / e2SquareLength > squaredMaxHeightFactor || e2SquareLength / e0SquareLength > squaredMaxHeightFactor )
|
|
{
|
|
double e1SquareLength = ( c2 - c1 ).lengthSquared();
|
|
double e3SquareLength = ( c0 - c3 ).lengthSquared();
|
|
|
|
if ( e1SquareLength / e3SquareLength > squaredMaxHeightFactor || e3SquareLength / e1SquareLength > squaredMaxHeightFactor )
|
|
{
|
|
return true;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool RigCell::isCollapsedCell( double nodeNearTolerance ) const
|
|
{
|
|
const std::vector<cvf::Vec3d>& nodes = m_hostGrid->mainGrid()->nodes();
|
|
|
|
cvf::ubyte faceVertexIndices[4];
|
|
cvf::ubyte oppFaceVertexIndices[4];
|
|
|
|
int face;
|
|
for ( face = 0; face < 6; face += 2 )
|
|
{
|
|
cvf::StructGridInterface::cellFaceVertexIndices( static_cast<cvf::StructGridInterface::FaceType>( face ), faceVertexIndices );
|
|
cvf::StructGridInterface::cellFaceVertexIndices( cvf::StructGridInterface::oppositeFace(
|
|
static_cast<cvf::StructGridInterface::FaceType>( face ) ),
|
|
oppFaceVertexIndices );
|
|
|
|
cvf::Vec3d c0 = nodes[m_cornerIndices[faceVertexIndices[0]]];
|
|
cvf::Vec3d c1 = nodes[m_cornerIndices[faceVertexIndices[1]]];
|
|
cvf::Vec3d c2 = nodes[m_cornerIndices[faceVertexIndices[2]]];
|
|
cvf::Vec3d c3 = nodes[m_cornerIndices[faceVertexIndices[3]]];
|
|
|
|
cvf::Vec3d oc0 = nodes[m_cornerIndices[oppFaceVertexIndices[0]]];
|
|
cvf::Vec3d oc1 = nodes[m_cornerIndices[oppFaceVertexIndices[1]]];
|
|
cvf::Vec3d oc2 = nodes[m_cornerIndices[oppFaceVertexIndices[2]]];
|
|
cvf::Vec3d oc3 = nodes[m_cornerIndices[oppFaceVertexIndices[3]]];
|
|
|
|
int zeroLengthEdgeCount = 0;
|
|
if ( isNear( c0, oc0, nodeNearTolerance ) )
|
|
{
|
|
++zeroLengthEdgeCount;
|
|
}
|
|
if ( isNear( c1, oc3, nodeNearTolerance ) )
|
|
{
|
|
++zeroLengthEdgeCount;
|
|
}
|
|
if ( isNear( c2, oc2, nodeNearTolerance ) )
|
|
{
|
|
++zeroLengthEdgeCount;
|
|
}
|
|
if ( isNear( c3, oc1, nodeNearTolerance ) )
|
|
{
|
|
++zeroLengthEdgeCount;
|
|
}
|
|
|
|
if ( zeroLengthEdgeCount >= 4 )
|
|
{
|
|
return true;
|
|
}
|
|
}
|
|
|
|
return false;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
cvf::Vec3d RigCell::faceCenter( cvf::StructGridInterface::FaceType face ) const
|
|
{
|
|
cvf::Vec3d avg( cvf::Vec3d::ZERO );
|
|
|
|
cvf::ubyte faceVertexIndices[4];
|
|
cvf::StructGridInterface::cellFaceVertexIndices( face, faceVertexIndices );
|
|
|
|
const std::vector<cvf::Vec3d>& nodeCoords = m_hostGrid->mainGrid()->nodes();
|
|
|
|
size_t i;
|
|
for ( i = 0; i < 4; i++ )
|
|
{
|
|
avg += nodeCoords[m_cornerIndices[faceVertexIndices[i]]];
|
|
}
|
|
|
|
avg /= 4.0;
|
|
|
|
return avg;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// Returns an area vector for the cell face. The direction is the face normal, and the length is
|
|
/// equal to the face area (projected to the plane represented by the diagonal in case of warp)
|
|
/// The components of this area vector are equal to the area of the face projection onto
|
|
/// the corresponding plane.
|
|
/// See http://geomalgorithms.com/a01-_area.html
|
|
//--------------------------------------------------------------------------------------------------
|
|
cvf::Vec3d RigCell::faceNormalWithAreaLength( cvf::StructGridInterface::FaceType face ) const
|
|
{
|
|
cvf::ubyte faceVertexIndices[4];
|
|
cvf::StructGridInterface::cellFaceVertexIndices( face, faceVertexIndices );
|
|
const std::vector<cvf::Vec3d>& nodeCoords = m_hostGrid->mainGrid()->nodes();
|
|
|
|
return 0.5 * ( nodeCoords[m_cornerIndices[faceVertexIndices[2]]] - nodeCoords[m_cornerIndices[faceVertexIndices[0]]] ) ^
|
|
( nodeCoords[m_cornerIndices[faceVertexIndices[3]]] - nodeCoords[m_cornerIndices[faceVertexIndices[1]]] );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigCell::volume() const
|
|
{
|
|
if ( m_isInvalid ) return 0.0;
|
|
const std::vector<cvf::Vec3d>& nodeCoords = m_hostGrid->mainGrid()->nodes();
|
|
|
|
std::array<cvf::Vec3d, 8> hexCorners;
|
|
for ( size_t i = 0; i < 8; ++i )
|
|
{
|
|
hexCorners[i] = nodeCoords.at( m_cornerIndices[i] );
|
|
}
|
|
return RigCellGeometryTools::calculateCellVolume( hexCorners );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// Find the intersection between the cell and the ray. The point closest to the ray origin is returned
|
|
/// in \a intersectionPoint, while the return value is the total number of intersections with the 24 triangles
|
|
/// the cell is interpreted as.
|
|
/// If no intersection is found, the intersection point is untouched.
|
|
//--------------------------------------------------------------------------------------------------
|
|
int RigCell::firstIntersectionPoint( const cvf::Ray& ray, cvf::Vec3d* intersectionPoint ) const
|
|
{
|
|
CVF_ASSERT( intersectionPoint != nullptr );
|
|
|
|
cvf::ubyte faceVertexIndices[4];
|
|
int face;
|
|
const std::vector<cvf::Vec3d>& nodes = m_hostGrid->mainGrid()->nodes();
|
|
|
|
cvf::Vec3d firstIntersection( cvf::Vec3d::ZERO );
|
|
double minLsq = HUGE_VAL;
|
|
int intersectionCount = 0;
|
|
|
|
for ( face = 0; face < 6; ++face )
|
|
{
|
|
cvf::StructGridInterface::cellFaceVertexIndices( static_cast<cvf::StructGridInterface::FaceType>( face ), faceVertexIndices );
|
|
cvf::Vec3d intersection;
|
|
cvf::Vec3d faceCenter = this->faceCenter( static_cast<cvf::StructGridInterface::FaceType>( face ) );
|
|
|
|
for ( size_t i = 0; i < 4; ++i )
|
|
{
|
|
size_t next = i < 3 ? i + 1 : 0;
|
|
if ( ray.triangleIntersect( nodes[m_cornerIndices[faceVertexIndices[i]]],
|
|
nodes[m_cornerIndices[faceVertexIndices[next]]],
|
|
faceCenter,
|
|
&intersection ) )
|
|
{
|
|
intersectionCount++;
|
|
double lsq = ( intersection - ray.origin() ).lengthSquared();
|
|
if ( lsq < minLsq )
|
|
{
|
|
firstIntersection = intersection;
|
|
minLsq = lsq;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
if ( intersectionCount > 0 )
|
|
{
|
|
*intersectionPoint = firstIntersection;
|
|
}
|
|
|
|
return intersectionCount;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RigCell::faceIndices( cvf::StructGridInterface::FaceType face, std::array<size_t, 4>* indices ) const
|
|
{
|
|
cvf::ubyte faceVertexIndices[4];
|
|
cvf::StructGridInterface::cellFaceVertexIndices( face, faceVertexIndices );
|
|
|
|
( *indices )[0] = m_cornerIndices[faceVertexIndices[0]];
|
|
( *indices )[1] = m_cornerIndices[faceVertexIndices[1]];
|
|
( *indices )[2] = m_cornerIndices[faceVertexIndices[2]];
|
|
( *indices )[3] = m_cornerIndices[faceVertexIndices[3]];
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
cvf::BoundingBox RigCell::boundingBox() const
|
|
{
|
|
cvf::BoundingBox bb;
|
|
std::array<cvf::Vec3d, 8> hexCorners;
|
|
|
|
if ( m_hostGrid )
|
|
{
|
|
m_hostGrid->cellCornerVertices( m_gridLocalCellIndex, hexCorners.data() );
|
|
for ( const auto& corner : hexCorners )
|
|
{
|
|
bb.add( corner );
|
|
}
|
|
}
|
|
return bb;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// Return the main grid neighbor cell of the given face
|
|
//--------------------------------------------------------------------------------------------------
|
|
RigCell RigCell::neighborCell( cvf::StructGridInterface::FaceType face ) const
|
|
{
|
|
size_t i, j, k;
|
|
|
|
m_hostGrid->ijkFromCellIndexUnguarded( mainGridCellIndex(), &i, &j, &k );
|
|
|
|
size_t neighborIdx;
|
|
if ( m_hostGrid->cellIJKNeighbor( i, j, k, face, &neighborIdx ) )
|
|
{
|
|
return m_hostGrid->cell( neighborIdx );
|
|
}
|
|
|
|
RigCell retcell;
|
|
retcell.setInvalid( true );
|
|
return retcell;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// Find and return the main grid cell index of all up to 26 neighbor cells to the given cell
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<size_t> RigCell::allNeighborMainGridCellIndexes() const
|
|
{
|
|
std::vector<size_t> neighbors;
|
|
|
|
size_t ni, nj, nk;
|
|
|
|
m_hostGrid->ijkFromCellIndexUnguarded( mainGridCellIndex(), &ni, &nj, &nk );
|
|
|
|
for ( size_t i = ni - 1; i <= ni + 1; i++ )
|
|
{
|
|
for ( size_t j = nj - 1; j <= nj + 1; j++ )
|
|
{
|
|
for ( size_t k = nk - 1; k <= nk + 1; k++ )
|
|
{
|
|
if ( m_hostGrid->isCellValid( i, j, k ) )
|
|
{
|
|
size_t cellIndex = m_hostGrid->cellIndexFromIJK( i, j, k );
|
|
if ( ( ni == i ) && ( nj == j ) && ( nk == k ) ) continue;
|
|
neighbors.push_back( cellIndex );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
return neighbors;
|
|
}
|