Files
ResInsight/ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartGrid.cpp
T
Magne Sjaastad 56e98ba1c4 #14476 ApplicationLibCode: Replace CVF_ASSERT with CAF_ASSERT
Migrate all assert macros in ApplicationLibCode to CAF_ASSERT and remove every
use of cvfAssert.h.

CVF_ASSERT is replaced one to one. CVF_TIGHT_ASSERT is also replaced by
CAF_ASSERT, which is semantically exact: CVF_ENABLE_TIGHT_ASSERTS is 1 only
under _DEBUG, and that is what CAF_ASSERT now does. The two CVF_FAIL_MSG sites
become CAF_ASSERT( false && "message" ), preserving the message with the idiom
already used elsewhere in the code base.

Counts before and after: CVF_ASSERT 1044 to 0, CVF_TIGHT_ASSERT 66 to 0,
CVF_FAIL_MSG 2 to 0, cvfAssert.h references 154 to 0.

Include handling: files that included cvfAssert.h directly now include
cafAssert.h instead, includes left dead by the migration are removed, and files
that were relying on cvfAssert.h transitively get an explicit cafAssert.h. Files
that reach cafAssert.h through another caf header are left unchanged; a missing
include here is a compile error, not a silently disabled assert.

ResultStatisticsCache links only LibCore and therefore had no path to
cafAssert.h. Add the cafPdmCore directory as a private include path rather than
linking the library, since cafAssert.h is header only.

Note that this stops these asserts from firing in Release and RelWithDebInfo,
where CVF_ASSERT was previously active.
2026-08-07 15:05:57 +02:00

535 lines
21 KiB
C++

/////////////////////////////////////////////////////////////////////////////////
//
// 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 "RigFemPartGrid.h"
#include "RigFemPart.h"
#include "cafAssert.h"
#include <array>
#include <climits>
#include <cmath>
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigFemPartGrid::RigFemPartGrid()
: m_femPart( nullptr )
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigFemPartGrid::~RigFemPartGrid()
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigFemPartGrid::setFemPart( const RigFemPart* femPart, bool invertIJK )
{
m_femPart = femPart;
generateStructGridData( invertIJK );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigFemPartGrid::generateStructGridData( bool invertIJK )
{
//[X] 1. Calculate neighbors for each element
//[X] record the ones with 3 or fewer neighbors as possible grid corners
//[X] 2. Loop over the possible corner cells,
//[X] find the one that corresponds to IJK = 000
//[X] by finding the one closest to origo // Does not work
//[X] by Determining what surfs correspond to NEG IJK surfaces in that element,
// and that none of those faces have a neighbor
//[X] 4. Assign IJK = 000 to that element
//[X] Store IJK in elm idx array
//[X] 5. Loop along POS I surfaces increment I for each element and assign IJK
//[X] when at end, go to POS J neighbor, increment J, repeat above.
//[X] etc for POS Z
//[X] Find max IJK as you go,
//[ ] also assert that there are no NEG I/NEG J/NEG Z neighbors when starting on a new row
//[ ] (Need to find min, and offset IJK values if there exists such)
//[ ] 6. If IJK to elm idx is needed, allocate "grid" with maxI,maxJ,maxZ values
//[ ] Loop over elms, assign elmIdx to IJK address in grid
int elmIdxForIJK_000 = findElmIdxForIJK000();
if ( elmIdxForIJK_000 == -1 ) return;
// Find the IJK faces based on the corner cell
cvf::Vec3i ijkMainFaceIndices = findMainIJKFaces( elmIdxForIJK_000 );
// assign ijk to cells
{
m_ijkPrElement.resize( m_femPart->elementCount(), cvf::Vec3i( -1, -1, -1 ) );
int posIFaceIdx = ijkMainFaceIndices[0];
int posJFaceIdx = ijkMainFaceIndices[1];
int posKFaceIdx = ijkMainFaceIndices[2];
m_elementIJKCounts = cvf::Vec3st( 0, 0, 0 );
int elmIdxInK = elmIdxForIJK_000;
int kCoord = 0;
while ( true )
{
int elmIdxInJ = elmIdxInK;
cvf::Vec3f startElmInKNormalJ = m_femPart->faceNormal( elmIdxInJ, posJFaceIdx );
cvf::Vec3f startElmInKNormalI = m_femPart->faceNormal( elmIdxInJ, posIFaceIdx );
int jCoord = 0;
while ( true )
{
int elmIdxInI = elmIdxInJ;
cvf::Vec3f startElmInJNormalI = m_femPart->faceNormal( elmIdxInI, posIFaceIdx );
int iCoord = 0;
while ( true )
{
CAF_ASSERT( elmIdxInI >= 0 && size_t( elmIdxInI ) < m_ijkPrElement.size() );
// Assign ijk coordinate
m_ijkPrElement[elmIdxInI] = cvf::Vec3i( iCoord, jCoord, kCoord );
++iCoord;
// Find neighbor and exit if at end
int neighborElmIdx = m_femPart->elementNeighbor( elmIdxInI, posIFaceIdx );
if ( neighborElmIdx == -1 ) break;
// Find the continuing face in the neighbor element (opposite of the neighbor face)
int neighborNegFaceIdx = m_femPart->neighborFace( elmIdxInI, posIFaceIdx );
RigElementType eType = m_femPart->elementType( neighborElmIdx );
posIFaceIdx = RigFemTypes::oppositeFace( eType, neighborNegFaceIdx );
// Step to neighbor
elmIdxInI = neighborElmIdx;
}
// Scoped to show that nothing bleeds further to K-loop
{
if ( iCoord > static_cast<int>( m_elementIJKCounts[0] ) ) m_elementIJKCounts[0] = iCoord;
++jCoord;
// Find neighbor and exit if at end
int neighborElmIdx = m_femPart->elementNeighbor( elmIdxInJ, posJFaceIdx );
if ( neighborElmIdx == -1 ) break;
// Find the continuing face in the neighbor element (opposite of the neighbor face)
int neighborNegFaceIdx = m_femPart->neighborFace( elmIdxInJ, posJFaceIdx );
RigElementType eType = m_femPart->elementType( neighborElmIdx );
posJFaceIdx = RigFemTypes::oppositeFace( eType, neighborNegFaceIdx );
// Now where is posIFace of the new J cell ?
posIFaceIdx = perpendicularFaceInDirection( startElmInJNormalI, neighborNegFaceIdx, neighborElmIdx );
// Step to neighbor
elmIdxInJ = neighborElmIdx;
}
}
{
if ( jCoord > static_cast<int>( m_elementIJKCounts[1] ) ) m_elementIJKCounts[1] = jCoord;
++kCoord;
// Find neighbor and exit if at end
int neighborElmIdx = m_femPart->elementNeighbor( elmIdxInK, posKFaceIdx );
if ( neighborElmIdx == -1 ) break;
// Find the continuing face in the neighbor element (opposite of the neighbor face)
int neighborNegFaceIdx = m_femPart->neighborFace( elmIdxInK, posKFaceIdx );
RigElementType eType = m_femPart->elementType( neighborElmIdx );
posKFaceIdx = RigFemTypes::oppositeFace( eType, neighborNegFaceIdx );
// Now where is posJFace of the new K cell ?
posJFaceIdx = perpendicularFaceInDirection( startElmInKNormalJ, neighborNegFaceIdx, neighborElmIdx );
posIFaceIdx = perpendicularFaceInDirection( startElmInKNormalI, neighborNegFaceIdx, neighborElmIdx );
// Step to neighbor
elmIdxInK = neighborElmIdx;
}
}
if ( kCoord > static_cast<int>( m_elementIJKCounts[2] ) ) m_elementIJKCounts[2] = kCoord;
}
m_elmIdxPrIJK.resize( m_elementIJKCounts[0], m_elementIJKCounts[1], m_elementIJKCounts[2] );
if ( invertIJK )
{
size_t iCount = m_elementIJKCounts[0];
size_t jCount = m_elementIJKCounts[1];
size_t kCount = m_elementIJKCounts[2];
for ( auto& ijk : m_ijkPrElement )
{
if ( ijk[0] >= 0 ) ijk[0] = static_cast<int>( iCount ) - 1 - ijk[0];
if ( ijk[1] >= 0 ) ijk[1] = static_cast<int>( jCount ) - 1 - ijk[1];
if ( ijk[2] >= 0 ) ijk[2] = static_cast<int>( kCount ) - 1 - ijk[2];
}
}
for ( int elmIdx = 0; elmIdx < m_femPart->elementCount(); ++elmIdx )
{
size_t i, j, k;
bool validIndex = ijkFromCellIndex( elmIdx, &i, &j, &k );
if ( validIndex )
{
m_elmIdxPrIJK.at( i, j, k ) = elmIdx;
}
}
// IJK bounding box
m_reservoirIJKBoundingBox.first = cvf::Vec3st( INT_MAX, INT_MAX, INT_MAX );
m_reservoirIJKBoundingBox.second = cvf::Vec3st( 0, 0, 0 );
cvf::Vec3st& min = m_reservoirIJKBoundingBox.first;
cvf::Vec3st& max = m_reservoirIJKBoundingBox.second;
for ( int elmIdx = 0; elmIdx < m_femPart->elementCount(); ++elmIdx )
{
RigElementType elementType = m_femPart->elementType( elmIdx );
size_t i, j, k;
bool validIndex = ijkFromCellIndex( elmIdx, &i, &j, &k );
if ( ( elementType == RigElementType::HEX8P ) && validIndex )
{
if ( i < min.x() ) min.x() = i;
if ( j < min.y() ) min.y() = j;
if ( k < min.z() ) min.z() = k;
if ( i > max.x() ) max.x() = i;
if ( j > max.y() ) max.y() = j;
if ( k > max.z() ) max.z() = k;
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
int RigFemPartGrid::findElmIdxForIJK000()
{
const std::vector<int>& possibleGridCorners = m_femPart->possibleGridCornerElements();
size_t possibleCornerCount = possibleGridCorners.size();
for ( size_t pcIdx = 0; pcIdx < possibleCornerCount; ++pcIdx )
{
int elmIdx = possibleGridCorners[pcIdx];
cvf::Vec3i ijkMainFaceIndices = findMainIJKFaces( elmIdx );
if ( m_femPart->elementNeighbor( elmIdx, ijkMainFaceIndices[0] ) != -1 &&
m_femPart->elementNeighbor( elmIdx, ijkMainFaceIndices[1] ) != -1 &&
m_femPart->elementNeighbor( elmIdx, ijkMainFaceIndices[2] ) != -1 )
{
return elmIdx;
}
}
return -1;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3i RigFemPartGrid::findMainIJKFaces( int elementIndex ) const
{
cvf::Vec3i ijkMainFaceIndices = cvf::Vec3i( -1, -1, -1 );
RigElementType eType = m_femPart->elementType( elementIndex );
int faceCount = RigFemTypes::elementFaceCount( eType );
std::vector<cvf::Vec3f> normals( faceCount );
for ( int faceIdx = 0; faceIdx < faceCount; ++faceIdx )
{
normals[faceIdx] = m_femPart->faceNormal( elementIndex, faceIdx );
}
// Record three independent main direction vectors for the element, and what face they are created from
cvf::Vec3f mainElmDirections[3];
int mainElmDirOriginFaces[3];
if ( RigFemTypes::is8NodeElement( eType ) )
{
mainElmDirections[0] = normals[0] - normals[1]; // To get a better "average" direction vector
mainElmDirections[1] = normals[2] - normals[3];
mainElmDirections[2] = normals[4] - normals[5];
mainElmDirOriginFaces[0] = 0;
mainElmDirOriginFaces[1] = 2;
mainElmDirOriginFaces[2] = 4;
}
else
{
mainElmDirections[0] = cvf::Vec3f::ZERO;
mainElmDirections[1] = cvf::Vec3f::ZERO;
mainElmDirections[2] = cvf::Vec3f::ZERO;
CAF_ASSERT( false );
}
mainElmDirections[0].normalize();
mainElmDirections[1].normalize();
mainElmDirections[2].normalize();
// Match the element main directions with best XYZ match (IJK respectively)
// Find the mainElmDirection with the largest component starting with Z
// and use that for the corresponding IJK direction.
// Find the Z (for K) first. Then select among the other two the Y (for J),
// and select the remaining for I
int mainElmDirectionIdxForIJK[3] = { -1, -1, -1 };
for ( int cIdx = 2; cIdx >= 0; --cIdx ) // Check Z first as it is more important
{
double maxAbsComp = -1.0;
int usedDir1 = -1;
int usedDir2 = -1;
for ( int dIdx = 0; dIdx < 3; ++dIdx )
{
if ( dIdx == usedDir1 || dIdx == usedDir2 ) continue;
float absComp = fabs( mainElmDirections[dIdx][cIdx] );
if ( absComp > maxAbsComp )
{
maxAbsComp = absComp;
mainElmDirectionIdxForIJK[cIdx] = dIdx;
}
}
if ( usedDir1 == -1 )
usedDir1 = mainElmDirectionIdxForIJK[cIdx];
else
usedDir2 = mainElmDirectionIdxForIJK[cIdx];
}
// Assign the correct face based on the main direction
ijkMainFaceIndices[0] = ( mainElmDirections[mainElmDirectionIdxForIJK[0]] * cvf::Vec3f::X_AXIS > 0 )
? mainElmDirOriginFaces[mainElmDirectionIdxForIJK[0]]
: RigFemTypes::oppositeFace( eType, mainElmDirOriginFaces[mainElmDirectionIdxForIJK[0]] );
ijkMainFaceIndices[1] = ( mainElmDirections[mainElmDirectionIdxForIJK[1]] * cvf::Vec3f::Y_AXIS > 0 )
? mainElmDirOriginFaces[mainElmDirectionIdxForIJK[1]]
: RigFemTypes::oppositeFace( eType, mainElmDirOriginFaces[mainElmDirectionIdxForIJK[1]] );
ijkMainFaceIndices[2] = ( mainElmDirections[mainElmDirectionIdxForIJK[2]] * -cvf::Vec3f::Z_AXIS > 0 )
? mainElmDirOriginFaces[mainElmDirectionIdxForIJK[2]]
: RigFemTypes::oppositeFace( eType, mainElmDirOriginFaces[mainElmDirectionIdxForIJK[2]] );
return ijkMainFaceIndices;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::pair<cvf::Vec3st, cvf::Vec3st> RigFemPartGrid::reservoirIJKBoundingBox() const
{
return m_reservoirIJKBoundingBox;
}
//--------------------------------------------------------------------------------------------------
/// Find the face that is not perpFaceIdx or its opposite, and has normal closest to direction
//--------------------------------------------------------------------------------------------------
int RigFemPartGrid::perpendicularFaceInDirection( cvf::Vec3f direction, int perpFaceIdx, int elmIdx )
{
RigElementType eType = m_femPart->elementType( elmIdx );
int faceCount = RigFemTypes::elementFaceCount( eType );
int oppFace = RigFemTypes::oppositeFace( eType, perpFaceIdx );
double minDiffSqLength = HUGE_VAL;
cvf::Vec3f faceNormal;
direction.normalize();
int bestFace = -1;
for ( int faceIdx = 0; faceIdx < faceCount; ++faceIdx )
{
if ( faceIdx == perpFaceIdx || faceIdx == oppFace ) continue;
faceNormal = m_femPart->faceNormal( elmIdx, faceIdx );
faceNormal.normalize();
float diffSqLength = ( direction - faceNormal ).lengthSquared();
if ( diffSqLength < minDiffSqLength )
{
bestFace = faceIdx;
minDiffSqLength = diffSqLength;
}
}
return bestFace;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartGrid::isCellValid( size_t i, size_t j, size_t k ) const
{
CAF_ASSERT( false );
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigFemPartGrid::minCoordinate() const
{
CAF_ASSERT( false );
return cvf::Vec3d::ZERO;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigFemPartGrid::maxCoordinate() const
{
CAF_ASSERT( false );
return cvf::Vec3d::ZERO;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartGrid::cellIJKNeighbor( size_t i, size_t j, size_t k, FaceType face, size_t* neighborCellIndex ) const
{
CAF_ASSERT( false );
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::cellIndexFromIJK( size_t i, size_t j, size_t k ) const
{
return m_elmIdxPrIJK.at( i, j, k );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartGrid::ijkFromCellIndex( size_t cellIndex, size_t* i, size_t* j, size_t* k ) const
{
if ( cellIndex < m_ijkPrElement.size() )
{
int signed_i = m_ijkPrElement[cellIndex][0];
int signed_j = m_ijkPrElement[cellIndex][1];
int signed_k = m_ijkPrElement[cellIndex][2];
if ( signed_i >= 0 && signed_j >= 0 && signed_k >= 0 )
{
*i = signed_i;
*j = signed_j;
*k = signed_k;
return true;
}
}
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartGrid::cellIJKFromCoordinate( const cvf::Vec3d& coord, size_t* i, size_t* j, size_t* k ) const
{
CAF_ASSERT( false );
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::array<cvf::Vec3d, 8> RigFemPartGrid::cellCornerVertices( size_t cellIndex ) const
{
const std::vector<cvf::Vec3f>& nodeCoords = m_femPart->nodes().coordinates;
const int* cornerIndices = m_femPart->connectivities( cellIndex );
std::array<cvf::Vec3d, 8> vertices;
for ( size_t i = 0; i < 8; ++i )
{
vertices[i] = cvf::Vec3d( nodeCoords[cornerIndices[i]] );
}
return vertices;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigFemPartGrid::cellCentroid( size_t cellIndex ) const
{
std::array<cvf::Vec3d, 8> cellVertices = cellCornerVertices( cellIndex );
cvf::Vec3d centroid( 0.0, 0.0, 0.0 );
for ( int i = 0; i < 8; ++i )
{
centroid += cellVertices[i];
}
return centroid / 8.0;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::cellCountI() const
{
return m_elementIJKCounts[0];
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::cellCountJ() const
{
return m_elementIJKCounts[1];
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::cellCountK() const
{
return m_elementIJKCounts[2];
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigFemPartGrid::cellMinMaxCordinates( size_t cellIndex, cvf::Vec3d* minCoordinate, cvf::Vec3d* maxCoordinate ) const
{
CAF_ASSERT( false );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::gridPointIndexFromIJK( size_t i, size_t j, size_t k ) const
{
CAF_ASSERT( false );
return cvf::UNDEFINED_SIZE_T;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigFemPartGrid::gridPointCoordinate( size_t i, size_t j, size_t k ) const
{
CAF_ASSERT( false );
return cvf::Vec3d::ZERO;
}