diff --git a/ApplicationLibCode/ProjectDataModel/GeoMech/RimGeoMechContourMapProjection.cpp b/ApplicationLibCode/ProjectDataModel/GeoMech/RimGeoMechContourMapProjection.cpp index 5e8d3946b8..a338b84b15 100644 --- a/ApplicationLibCode/ProjectDataModel/GeoMech/RimGeoMechContourMapProjection.cpp +++ b/ApplicationLibCode/ProjectDataModel/GeoMech/RimGeoMechContourMapProjection.cpp @@ -23,6 +23,7 @@ #include "RiaWeightedMeanCalculator.h" #include "RigCellGeometryTools.h" +#include "RigContourMapGrid.h" #include "RigFemAddressDefines.h" #include "RigFemPart.h" #include "RigFemPartCollection.h" @@ -236,14 +237,15 @@ void RimGeoMechContourMapProjection::updateGridInformation() m_kLayers = m_femPartGrid->cellCountK(); m_femPart->ensureIntersectionSearchTreeIsBuilt(); - m_gridBoundingBox = geoMechCase->activeCellsBoundingBox(); + cvf::BoundingBox gridBoundingBox = geoMechCase->activeCellsBoundingBox(); + cvf::BoundingBox expandedBoundingBox; if ( m_limitToPorePressureRegions ) { auto [stepIdx, frameIdx] = view()->currentStepAndDataFrame(); - m_expandedBoundingBox = calculateExpandedPorBarBBox( stepIdx, frameIdx ); - if ( !m_expandedBoundingBox.isValid() ) + expandedBoundingBox = calculateExpandedPorBarBBox( stepIdx, frameIdx ); + if ( !expandedBoundingBox.isValid() ) { m_limitToPorePressureRegions = false; } @@ -251,26 +253,19 @@ void RimGeoMechContourMapProjection::updateGridInformation() if ( !m_limitToPorePressureRegions ) { - m_expandedBoundingBox = m_gridBoundingBox; + expandedBoundingBox = gridBoundingBox; } - cvf::Vec3d minExpandedPoint = m_expandedBoundingBox.min() - cvf::Vec3d( gridEdgeOffset(), gridEdgeOffset(), 0.0 ); - cvf::Vec3d maxExpandedPoint = m_expandedBoundingBox.max() + cvf::Vec3d( gridEdgeOffset(), gridEdgeOffset(), 0.0 ); + cvf::Vec3d minExpandedPoint = expandedBoundingBox.min() - cvf::Vec3d( gridEdgeOffset(), gridEdgeOffset(), 0.0 ); + cvf::Vec3d maxExpandedPoint = expandedBoundingBox.max() + cvf::Vec3d( gridEdgeOffset(), gridEdgeOffset(), 0.0 ); if ( m_limitToPorePressureRegions && !m_applyPPRegionLimitVertically ) { - minExpandedPoint.z() = m_gridBoundingBox.min().z(); - maxExpandedPoint.z() = m_gridBoundingBox.max().z(); + minExpandedPoint.z() = gridBoundingBox.min().z(); + maxExpandedPoint.z() = gridBoundingBox.max().z(); } - m_expandedBoundingBox = cvf::BoundingBox( minExpandedPoint, maxExpandedPoint ); + expandedBoundingBox = cvf::BoundingBox( minExpandedPoint, maxExpandedPoint ); - m_mapSize = calculateMapSize(); - - // Re-jig max point to be an exact multiple of cell size - cvf::Vec3d minPoint = m_expandedBoundingBox.min(); - cvf::Vec3d maxPoint = m_expandedBoundingBox.max(); - maxPoint.x() = minPoint.x() + m_mapSize.x() * sampleSpacing(); - maxPoint.y() = minPoint.y() + m_mapSize.y() * sampleSpacing(); - m_expandedBoundingBox = cvf::BoundingBox( minPoint, maxPoint ); + m_contourMapGrid = std::make_unique( gridBoundingBox, expandedBoundingBox, sampleSpacing() ); } //-------------------------------------------------------------------------------------------------- @@ -298,11 +293,11 @@ std::vector RimGeoMechContourMapProjection::getMapCellVisibility() cvf::BoundingBox validResBoundingBox; for ( size_t cellIndex = 0; cellIndex < cellResults.size(); ++cellIndex ) { - cvf::Vec2ui ij = ijFromCellIndex( cellIndex ); + cvf::Vec2ui ij = m_contourMapGrid->ijFromCellIndex( cellIndex ); if ( cellResults[cellIndex] != std::numeric_limits::infinity() ) { distanceImage[ij.x()][ij.y()] = 1u; - validResBoundingBox.add( cvf::Vec3d( cellCenterPosition( ij.x(), ij.y() ), 0.0 ) ); + validResBoundingBox.add( cvf::Vec3d( m_contourMapGrid->cellCenterPosition( ij.x(), ij.y() ), 0.0 ) ); } else { @@ -322,7 +317,7 @@ std::vector RimGeoMechContourMapProjection::getMapCellVisibility() { if ( !mapCellVisibility[cellIndex] ) { - cvf::Vec2ui ij = ijFromCellIndex( cellIndex ); + cvf::Vec2ui ij = m_contourMapGrid->ijFromCellIndex( cellIndex ); if ( distanceImage[ij.x()][ij.y()] < cellPadding * cellPadding ) { mapCellVisibility[cellIndex] = true; @@ -412,7 +407,7 @@ std::vector RimGeoMechContourMapProjection::generateResultsFromAddress( { if ( mapCellVisibility.empty() || mapCellVisibility[index] ) { - cvf::Vec2ui ij = ijFromCellIndex( index ); + cvf::Vec2ui ij = m_contourMapGrid->ijFromCellIndex( index ); aggregatedResults[index] = calculateValueInMapCell( ij.x(), ij.y(), resultValues ); } } diff --git a/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.cpp b/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.cpp index 793286a413..18e27b50c9 100644 --- a/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.cpp +++ b/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.cpp @@ -24,6 +24,7 @@ #include "RiaWeightedMeanCalculator.h" #include "RigCellGeometryTools.h" +#include "RigContourMapGrid.h" #include "RigHexIntersectionTools.h" #include "RimCase.h" @@ -75,7 +76,6 @@ CAF_PDM_ABSTRACT_SOURCE_INIT( RimContourMapProjection, "RimContourMapProjection" //-------------------------------------------------------------------------------------------------- RimContourMapProjection::RimContourMapProjection() : m_pickPoint( cvf::Vec2d::UNDEFINED ) - , m_mapSize( cvf::Vec2ui( 0u, 0u ) ) , m_currentResultTimestep( -1 ) , m_minResultAllTimeSteps( std::numeric_limits::infinity() ) , m_maxResultAllTimeSteps( -std::numeric_limits::infinity() ) @@ -166,8 +166,8 @@ std::vector RimContourMapProjection::generatePickPointPolygon() { #ifndef NDEBUG cvf::Vec2d cellDiagonal( sampleSpacing() * 0.5, sampleSpacing() * 0.5 ); - cvf::Vec2ui pickedCell = ijFromLocalPos( m_pickPoint ); - cvf::Vec2d cellCenter = cellCenterPosition( pickedCell.x(), pickedCell.y() ); + cvf::Vec2ui pickedCell = m_contourMapGrid->ijFromLocalPos( m_pickPoint ); + cvf::Vec2d cellCenter = m_contourMapGrid->cellCenterPosition( pickedCell.x(), pickedCell.y() ); cvf::Vec2d cellCorner = cellCenter - cellDiagonal; points.push_back( cvf::Vec3d( cellCorner, 0.0 ) ); points.push_back( cvf::Vec3d( cellCorner + cvf::Vec2d( sampleSpacing(), 0.0 ), 0.0 ) ); @@ -326,7 +326,7 @@ double RimContourMapProjection::sumAllValues() const //-------------------------------------------------------------------------------------------------- cvf::Vec2ui RimContourMapProjection::numberOfElementsIJ() const { - return m_mapSize; + return m_contourMapGrid->numberOfElementsIJ(); } //-------------------------------------------------------------------------------------------------- @@ -334,10 +334,7 @@ cvf::Vec2ui RimContourMapProjection::numberOfElementsIJ() const //-------------------------------------------------------------------------------------------------- cvf::Vec2ui RimContourMapProjection::numberOfVerticesIJ() const { - cvf::Vec2ui mapSize = numberOfElementsIJ(); - mapSize.x() += 1u; - mapSize.y() += 1u; - return mapSize; + return m_contourMapGrid->numberOfVerticesIJ(); } //-------------------------------------------------------------------------------------------------- @@ -354,7 +351,7 @@ bool RimContourMapProjection::isColumnResult() const //-------------------------------------------------------------------------------------------------- double RimContourMapProjection::valueAtVertex( uint i, uint j ) const { - size_t index = vertexIndexFromIJ( i, j ); + size_t index = m_contourMapGrid->vertexIndexFromIJ( i, j ); if ( index < numberOfVertices() ) { return m_aggregatedVertexResults.at( index ); @@ -367,7 +364,7 @@ double RimContourMapProjection::valueAtVertex( uint i, uint j ) const //-------------------------------------------------------------------------------------------------- uint RimContourMapProjection::numberOfCells() const { - return m_mapSize.x() * m_mapSize.y(); + return m_contourMapGrid->numberOfCells(); } //-------------------------------------------------------------------------------------------------- @@ -378,7 +375,7 @@ uint RimContourMapProjection::numberOfValidCells() const uint validCount = 0u; for ( uint i = 0; i < numberOfCells(); ++i ) { - cvf::Vec2ui ij = ijFromCellIndex( i ); + cvf::Vec2ui ij = m_contourMapGrid->ijFromCellIndex( i ); if ( hasResultInCell( ij.x(), ij.y() ) ) { validCount++; @@ -404,9 +401,10 @@ bool RimContourMapProjection::checkForMapIntersection( const cvf::Vec3d& domainP CVF_TIGHT_ASSERT( contourMapPoint ); CVF_TIGHT_ASSERT( valueAtPoint ); - cvf::Vec3d mapPos3d = domainPoint3d - m_expandedBoundingBox.min(); - cvf::Vec2d mapPos2d( mapPos3d.x(), mapPos3d.y() ); - cvf::Vec2d gridorigin( m_expandedBoundingBox.min().x(), m_expandedBoundingBox.min().y() ); + const cvf::Vec3d& minPoint = m_contourMapGrid->expandedBoundingBox().min(); + cvf::Vec3d mapPos3d = domainPoint3d - minPoint; + cvf::Vec2d mapPos2d( mapPos3d.x(), mapPos3d.y() ); + cvf::Vec2d gridorigin( minPoint.x(), minPoint.y() ); double value = interpolateValue( mapPos2d ); if ( value != std::numeric_limits::infinity() ) @@ -424,7 +422,7 @@ bool RimContourMapProjection::checkForMapIntersection( const cvf::Vec3d& domainP //-------------------------------------------------------------------------------------------------- void RimContourMapProjection::setPickPoint( cvf::Vec2d globalPickPoint ) { - m_pickPoint = globalPickPoint - origin2d(); + m_pickPoint = globalPickPoint - m_contourMapGrid->origin2d(); } //-------------------------------------------------------------------------------------------------- @@ -432,7 +430,7 @@ void RimContourMapProjection::setPickPoint( cvf::Vec2d globalPickPoint ) //-------------------------------------------------------------------------------------------------- cvf::Vec3d RimContourMapProjection::origin3d() const { - return m_expandedBoundingBox.min(); + return m_contourMapGrid->expandedBoundingBox().min(); } //-------------------------------------------------------------------------------------------------- @@ -774,9 +772,9 @@ std::vector>> RimContourMapProjection::gen #pragma omp parallel for for ( int index = 0; index < nCells; ++index ) { - cvf::Vec2ui ij = ijFromCellIndex( index ); + cvf::Vec2ui ij = m_contourMapGrid->ijFromCellIndex( index ); - cvf::Vec2d globalPos = cellCenterPosition( ij.x(), ij.y() ) + origin2d(); + cvf::Vec2d globalPos = m_contourMapGrid->cellCenterPosition( ij.x(), ij.y() ) + m_contourMapGrid->origin2d(); projected3dGridIndices[index] = cellRayIntersectionAndResults( globalPos, weightingResultValues ); } } @@ -785,9 +783,9 @@ std::vector>> RimContourMapProjection::gen #pragma omp parallel for for ( int index = 0; index < nCells; ++index ) { - cvf::Vec2ui ij = ijFromCellIndex( index ); + cvf::Vec2ui ij = m_contourMapGrid->ijFromCellIndex( index ); - cvf::Vec2d globalPos = cellCenterPosition( ij.x(), ij.y() ) + origin2d(); + cvf::Vec2d globalPos = m_contourMapGrid->cellCenterPosition( ij.x(), ij.y() ) + m_contourMapGrid->origin2d(); projected3dGridIndices[index] = cellOverlapVolumesAndResults( globalPos, weightingResultValues ); } } @@ -808,7 +806,7 @@ void RimContourMapProjection::generateVertexResults() #pragma omp parallel for for ( int index = 0; index < static_cast( nVertices ); ++index ) { - cvf::Vec2ui ij = ijFromVertexIndex( index ); + cvf::Vec2ui ij = m_contourMapGrid->ijFromVertexIndex( index ); m_aggregatedVertexResults[index] = calculateValueAtVertex( ij.x(), ij.y() ); } } @@ -818,7 +816,7 @@ void RimContourMapProjection::generateVertexResults() //-------------------------------------------------------------------------------------------------- void RimContourMapProjection::generateTrianglesWithVertexValues() { - std::vector vertices = generateVertices(); + std::vector vertices = m_contourMapGrid->generateVertices(); cvf::Vec2ui patchSize = numberOfVerticesIJ(); cvf::ref faceList = new cvf::UIntArray; @@ -1045,29 +1043,6 @@ void RimContourMapProjection::generateTrianglesWithVertexValues() m_trianglesWithVertexValues = finalTriangles; } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -std::vector RimContourMapProjection::generateVertices() const -{ - size_t nVertices = numberOfVertices(); - std::vector vertices( nVertices, cvf::Vec3d::ZERO ); - -#pragma omp parallel for - for ( int index = 0; index < static_cast( nVertices ); ++index ) - { - cvf::Vec2ui ij = ijFromVertexIndex( index ); - cvf::Vec2d mapPos = cellCenterPosition( ij.x(), ij.y() ); - // Shift away from sample point to vertex - mapPos.x() -= sampleSpacing() * 0.5; - mapPos.y() -= sampleSpacing() * 0.5; - - cvf::Vec3d vertexPos( mapPos, 0.0 ); - vertices[index] = vertexPos; - } - return vertices; -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -1301,18 +1276,19 @@ double RimContourMapProjection::sumTriangleAreas( const std::vector& std::vector RimContourMapProjection::cellOverlapVolumesAndResults( const cvf::Vec2d& globalPos2d, const std::vector& weightingResultValues ) const { - cvf::Vec3d top2dElementCentroid( globalPos2d, m_expandedBoundingBox.max().z() ); - cvf::Vec3d bottom2dElementCentroid( globalPos2d, m_expandedBoundingBox.min().z() ); - cvf::Vec3d planarDiagonalVector( 0.5 * sampleSpacing(), 0.5 * sampleSpacing(), 0.0 ); - cvf::Vec3d topNECorner = top2dElementCentroid + planarDiagonalVector; - cvf::Vec3d bottomSWCorner = bottom2dElementCentroid - planarDiagonalVector; + const cvf::BoundingBox& expandedBoundingBox = m_contourMapGrid->expandedBoundingBox(); + cvf::Vec3d top2dElementCentroid( globalPos2d, expandedBoundingBox.max().z() ); + cvf::Vec3d bottom2dElementCentroid( globalPos2d, expandedBoundingBox.min().z() ); + cvf::Vec3d planarDiagonalVector( 0.5 * sampleSpacing(), 0.5 * sampleSpacing(), 0.0 ); + cvf::Vec3d topNECorner = top2dElementCentroid + planarDiagonalVector; + cvf::Vec3d bottomSWCorner = bottom2dElementCentroid - planarDiagonalVector; cvf::BoundingBox bbox2dElement( bottomSWCorner, topNECorner ); std::vector> matchingVisibleCellsAndWeight; // Bounding box has been expanded, so 2d element may be outside actual 3d grid - if ( !bbox2dElement.intersects( m_gridBoundingBox ) ) + if ( !bbox2dElement.intersects( m_contourMapGrid->originalBoundingBox() ) ) { return matchingVisibleCellsAndWeight; } @@ -1362,11 +1338,13 @@ std::vector { std::vector> matchingVisibleCellsAndWeight; - cvf::Vec3d highestPoint( globalPos2d, m_expandedBoundingBox.max().z() ); - cvf::Vec3d lowestPoint( globalPos2d, m_expandedBoundingBox.min().z() ); + const cvf::BoundingBox& expandedBoundingBox = m_contourMapGrid->expandedBoundingBox(); + + cvf::Vec3d highestPoint( globalPos2d, expandedBoundingBox.max().z() ); + cvf::Vec3d lowestPoint( globalPos2d, expandedBoundingBox.min().z() ); // Bounding box has been expanded, so ray may be outside actual 3d grid - if ( !m_gridBoundingBox.contains( highestPoint ) ) + if ( !m_contourMapGrid->originalBoundingBox().contains( highestPoint ) ) { return matchingVisibleCellsAndWeight; } @@ -1441,8 +1419,8 @@ bool RimContourMapProjection::isStraightSummationResult( ResultAggregationEnum a //-------------------------------------------------------------------------------------------------- double RimContourMapProjection::interpolateValue( const cvf::Vec2d& gridPos2d ) const { - cvf::Vec2ui cellContainingPoint = ijFromLocalPos( gridPos2d ); - cvf::Vec2d cellCenter = cellCenterPosition( cellContainingPoint.x(), cellContainingPoint.y() ); + cvf::Vec2ui cellContainingPoint = m_contourMapGrid->ijFromLocalPos( gridPos2d ); + cvf::Vec2d cellCenter = m_contourMapGrid->cellCenterPosition( cellContainingPoint.x(), cellContainingPoint.y() ); std::array x; x[0] = cvf::Vec3d( cellCenter + cvf::Vec2d( -sampleSpacing() * 0.5, -sampleSpacing() * 0.5 ), 0.0 ); @@ -1494,7 +1472,7 @@ double RimContourMapProjection::interpolateValue( const cvf::Vec2d& gridPos2d ) //-------------------------------------------------------------------------------------------------- double RimContourMapProjection::valueInCell( uint i, uint j ) const { - size_t index = cellIndexFromIJ( i, j ); + size_t index = m_contourMapGrid->cellIndexFromIJ( i, j ); if ( index < numberOfCells() ) { return m_aggregatedResults.at( index ); @@ -1520,8 +1498,8 @@ double RimContourMapProjection::calculateValueAtVertex( uint vi, uint vj ) const if ( vi > 0u ) averageIs.push_back( vi - 1 ); if ( vj > 0u ) averageJs.push_back( vj - 1 ); - if ( vi < m_mapSize.x() ) averageIs.push_back( vi ); - if ( vj < m_mapSize.y() ) averageJs.push_back( vj ); + if ( vi < m_contourMapGrid->mapSize().x() ) averageIs.push_back( vi ); + if ( vj < m_contourMapGrid->mapSize().y() ) averageJs.push_back( vj ); RiaWeightedMeanCalculator calc; for ( uint j : averageJs ) @@ -1546,7 +1524,7 @@ double RimContourMapProjection::calculateValueAtVertex( uint vi, uint vj ) const //-------------------------------------------------------------------------------------------------- std::vector> RimContourMapProjection::cellsAtIJ( uint i, uint j ) const { - size_t cellIndex = cellIndexFromIJ( i, j ); + size_t cellIndex = m_contourMapGrid->cellIndexFromIJ( i, j ); if ( cellIndex < m_projected3dGridIndices.size() ) { return m_projected3dGridIndices[cellIndex]; @@ -1554,96 +1532,12 @@ std::vector> RimContourMapProjection::cellsAtIJ( uint return std::vector>(); } -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -size_t RimContourMapProjection::cellIndexFromIJ( uint i, uint j ) const -{ - CVF_ASSERT( i < m_mapSize.x() ); - CVF_ASSERT( j < m_mapSize.y() ); - - return i + j * m_mapSize.x(); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -size_t RimContourMapProjection::vertexIndexFromIJ( uint i, uint j ) const -{ - return i + j * ( m_mapSize.x() + 1 ); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -cvf::Vec2ui RimContourMapProjection::ijFromVertexIndex( size_t gridIndex ) const -{ - cvf::Vec2ui gridSize = numberOfVerticesIJ(); - - uint quotientX = static_cast( gridIndex ) / gridSize.x(); - uint remainderX = static_cast( gridIndex ) % gridSize.x(); - - return cvf::Vec2ui( remainderX, quotientX ); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -cvf::Vec2ui RimContourMapProjection::ijFromCellIndex( size_t cellIndex ) const -{ - CVF_TIGHT_ASSERT( cellIndex < numberOfCells() ); - - uint quotientX = static_cast( cellIndex ) / m_mapSize.x(); - uint remainderX = static_cast( cellIndex ) % m_mapSize.x(); - - return cvf::Vec2ui( remainderX, quotientX ); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -cvf::Vec2ui RimContourMapProjection::ijFromLocalPos( const cvf::Vec2d& localPos2d ) const -{ - uint i = localPos2d.x() / sampleSpacing(); - uint j = localPos2d.y() / sampleSpacing(); - return cvf::Vec2ui( i, j ); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -cvf::Vec2d RimContourMapProjection::cellCenterPosition( uint i, uint j ) const -{ - cvf::Vec3d gridExtent = m_expandedBoundingBox.extent(); - cvf::Vec2d cellCorner = cvf::Vec2d( ( i * gridExtent.x() ) / ( m_mapSize.x() ), ( j * gridExtent.y() ) / ( m_mapSize.y() ) ); - - return cellCorner + cvf::Vec2d( sampleSpacing() * 0.5, sampleSpacing() * 0.5 ); -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -cvf::Vec2d RimContourMapProjection::origin2d() const -{ - return cvf::Vec2d( m_expandedBoundingBox.min().x(), m_expandedBoundingBox.min().y() ); -} - //-------------------------------------------------------------------------------------------------- /// Vertex positions in local coordinates (add origin2d.x() for UTM x) //-------------------------------------------------------------------------------------------------- std::vector RimContourMapProjection::xVertexPositions() const { - double gridExtent = m_expandedBoundingBox.extent().x(); - - cvf::Vec2ui gridSize = numberOfVerticesIJ(); - std::vector positions; - positions.reserve( gridSize.x() ); - for ( uint i = 0; i < gridSize.x(); ++i ) - { - positions.push_back( ( i * gridExtent ) / ( gridSize.x() - 1 ) ); - } - - return positions; + return m_contourMapGrid->xVertexPositions(); } //-------------------------------------------------------------------------------------------------- @@ -1651,30 +1545,7 @@ std::vector RimContourMapProjection::xVertexPositions() const //-------------------------------------------------------------------------------------------------- std::vector RimContourMapProjection::yVertexPositions() const { - double gridExtent = m_expandedBoundingBox.extent().y(); - - cvf::Vec2ui gridSize = numberOfVerticesIJ(); - std::vector positions; - positions.reserve( gridSize.y() ); - for ( uint j = 0; j < gridSize.y(); ++j ) - { - positions.push_back( ( j * gridExtent ) / ( gridSize.y() - 1 ) ); - } - - return positions; -} - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -cvf::Vec2ui RimContourMapProjection::calculateMapSize() const -{ - cvf::Vec3d gridExtent = m_expandedBoundingBox.extent(); - - uint projectionSizeX = static_cast( std::ceil( gridExtent.x() / sampleSpacing() ) ); - uint projectionSizeY = static_cast( std::ceil( gridExtent.y() / sampleSpacing() ) ); - - return cvf::Vec2ui( projectionSizeX, projectionSizeY ); + return m_contourMapGrid->yVertexPositions(); } //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.h b/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.h index a0872066e5..c490291ece 100644 --- a/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.h +++ b/ApplicationLibCode/ProjectDataModel/RimContourMapProjection.h @@ -33,6 +33,7 @@ #include "cvfString.h" #include "cvfVector2.h" +class RigContourMapGrid; class RimGridView; class RimRegularLegendConfig; @@ -165,10 +166,9 @@ class RimContourMapProjection : public RimCheckableNamedObject bool mapCellVisibilityNeedsUpdating(); std::vector>> generateGridMapping(); - void generateVertexResults(); - void generateTrianglesWithVertexValues(); - std::vector generateVertices() const; - void generateContourPolygons(); + void generateVertexResults(); + void generateTrianglesWithVertexValues(); + void generateContourPolygons(); ContourPolygons createContourPolygonsFromLineSegments( caf::ContourLines::ListOfLineSegments& unorderedLineSegments, double contourValue ); void smoothContourPolygons( ContourPolygons* contourPolygons, bool favourExpansion ); void clipContourPolygons( ContourPolygons* contourPolygons, const ContourPolygons* clipBy ); @@ -191,15 +191,7 @@ class RimContourMapProjection : public RimCheckableNamedObject // Cell index and position conversion std::vector cellsAtIJ( uint i, uint j ) const; - size_t cellIndexFromIJ( uint i, uint j ) const; - size_t vertexIndexFromIJ( uint i, uint j ) const; - cvf::Vec2ui ijFromVertexIndex( size_t gridIndex ) const; - cvf::Vec2ui ijFromCellIndex( size_t mapIndex ) const; - cvf::Vec2ui ijFromLocalPos( const cvf::Vec2d& localPos2d ) const; - cvf::Vec2d cellCenterPosition( uint i, uint j ) const; - cvf::Vec2d origin2d() const; - - cvf::Vec2ui calculateMapSize() const; + double gridEdgeOffset() const; virtual void updateAfterResultGeneration( int timeStep ) = 0; @@ -225,9 +217,6 @@ class RimContourMapProjection : public RimCheckableNamedObject std::vector>> m_projected3dGridIndices; cvf::Vec2d m_pickPoint; - cvf::Vec2ui m_mapSize; - cvf::BoundingBox m_expandedBoundingBox; - cvf::BoundingBox m_gridBoundingBox; std::vector m_contourPolygons; std::vector m_contourLevelCumulativeAreas; std::vector m_trianglesWithVertexValues; @@ -236,4 +225,6 @@ class RimContourMapProjection : public RimCheckableNamedObject double m_minResultAllTimeSteps; double m_maxResultAllTimeSteps; + + std::unique_ptr m_contourMapGrid; }; diff --git a/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.cpp b/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.cpp index 5fd2902241..8b615f1cb5 100644 --- a/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.cpp +++ b/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.cpp @@ -26,6 +26,7 @@ #include "RigCaseCellResultsData.h" #include "RigCell.h" #include "RigCellGeometryTools.h" +#include "RigContourMapGrid.h" #include "RigEclipseCaseData.h" #include "RigEclipseResultAddress.h" #include "RigHexIntersectionTools.h" @@ -238,7 +239,7 @@ std::vector RimEclipseContourMapProjection::generateResults( int timeSte #pragma omp parallel for for ( int index = 0; index < static_cast( nCells ); ++index ) { - cvf::Vec2ui ij = ijFromCellIndex( index ); + cvf::Vec2ui ij = m_contourMapGrid->ijFromCellIndex( index ); aggregatedResults[index] = calculateValueInMapCell( ij.x(), ij.y(), gridResultValues ); } } @@ -343,19 +344,8 @@ void RimEclipseContourMapProjection::updateGridInformation() m_activeCellInfo = eclipseCase->eclipseCaseData()->activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL ); m_kLayers = m_mainGrid->cellCountK(); - m_gridBoundingBox = eclipseCase->activeCellsBoundingBox(); - cvf::Vec3d minExpandedPoint = m_gridBoundingBox.min() - cvf::Vec3d( gridEdgeOffset(), gridEdgeOffset(), 0.0 ); - cvf::Vec3d maxExpandedPoint = m_gridBoundingBox.max() + cvf::Vec3d( gridEdgeOffset(), gridEdgeOffset(), 0.0 ); - m_expandedBoundingBox = cvf::BoundingBox( minExpandedPoint, maxExpandedPoint ); - - m_mapSize = calculateMapSize(); - - // Re-jig max point to be an exact multiple of cell size - cvf::Vec3d minPoint = m_expandedBoundingBox.min(); - cvf::Vec3d maxPoint = m_expandedBoundingBox.max(); - maxPoint.x() = minPoint.x() + m_mapSize.x() * sampleSpacing(); - maxPoint.y() = minPoint.y() + m_mapSize.y() * sampleSpacing(); - m_expandedBoundingBox = cvf::BoundingBox( minPoint, maxPoint ); + cvf::BoundingBox gridBoundingBox = eclipseCase->activeCellsBoundingBox(); + m_contourMapGrid = std::make_unique( gridBoundingBox, sampleSpacing() ); } //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.h b/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.h index f3f1250a3b..4ec2333813 100644 --- a/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.h +++ b/ApplicationLibCode/ProjectDataModel/RimEclipseContourMapProjection.h @@ -32,6 +32,7 @@ class RigActiveCellInfo; class RigMainGrid; +class RigContourMapGrid; class RigResultAccessor; class RimEclipseContourMapView; class RimEclipseCase; diff --git a/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake b/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake index 574795281c..be2bd970b6 100644 --- a/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake +++ b/ApplicationLibCode/ReservoirDataModel/CMakeLists_files.cmake @@ -102,6 +102,7 @@ set(SOURCE_GROUP_HEADER_FILES ${CMAKE_CURRENT_LIST_DIR}/RigVfpTables.h ${CMAKE_CURRENT_LIST_DIR}/RigOsduWellLogData.h ${CMAKE_CURRENT_LIST_DIR}/RigWellTargetCandidatesGenerator.h + ${CMAKE_CURRENT_LIST_DIR}/RigContourMapGrid.h ) set(SOURCE_GROUP_SOURCE_FILES @@ -202,6 +203,7 @@ set(SOURCE_GROUP_SOURCE_FILES ${CMAKE_CURRENT_LIST_DIR}/RigVfpTables.cpp ${CMAKE_CURRENT_LIST_DIR}/RigOsduWellLogData.cpp ${CMAKE_CURRENT_LIST_DIR}/RigWellTargetCandidatesGenerator.cpp + ${CMAKE_CURRENT_LIST_DIR}/RigContourMapGrid.cpp ) list(APPEND CODE_HEADER_FILES ${SOURCE_GROUP_HEADER_FILES}) diff --git a/ApplicationLibCode/ReservoirDataModel/RigContourMapGrid.cpp b/ApplicationLibCode/ReservoirDataModel/RigContourMapGrid.cpp new file mode 100644 index 0000000000..c66502d478 --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/RigContourMapGrid.cpp @@ -0,0 +1,295 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2024- Equinor ASA +// +// 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 +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#include "RigContourMapGrid.h" + +#include + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigContourMapGrid::RigContourMapGrid( const cvf::BoundingBox& originalBoundingBox, double sampleSpacing ) + : m_sampleSpacing( sampleSpacing ) + , m_mapSize( cvf::Vec2ui( 0u, 0u ) ) + , m_originalBoundingBox( originalBoundingBox ) +{ + cvf::Vec3d minExpandedPoint = originalBoundingBox.min() - cvf::Vec3d( gridEdgeOffset(), gridEdgeOffset(), 0.0 ); + cvf::Vec3d maxExpandedPoint = originalBoundingBox.max() + cvf::Vec3d( gridEdgeOffset(), gridEdgeOffset(), 0.0 ); + m_expandedBoundingBox = cvf::BoundingBox( minExpandedPoint, maxExpandedPoint ); + + m_mapSize = calculateMapSize( m_expandedBoundingBox.extent(), sampleSpacing ); + + // Re-jig max point to be an exact multiple of cell size + m_expandedBoundingBox = makeMaxPointMultipleOfCellSize( m_expandedBoundingBox, m_mapSize, sampleSpacing ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigContourMapGrid::RigContourMapGrid( const cvf::BoundingBox& originalBoundingBox, const cvf::BoundingBox& expandedBoundingBox, double sampleSpacing ) + : m_sampleSpacing( sampleSpacing ) + , m_mapSize( cvf::Vec2ui( 0u, 0u ) ) + , m_originalBoundingBox( originalBoundingBox ) + , m_expandedBoundingBox( expandedBoundingBox ) +{ + m_mapSize = calculateMapSize( m_expandedBoundingBox.extent(), sampleSpacing ); + + // Re-jig max point to be an exact multiple of cell size + m_expandedBoundingBox = makeMaxPointMultipleOfCellSize( m_expandedBoundingBox, m_mapSize, sampleSpacing ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigContourMapGrid::sampleSpacing() const +{ + return m_sampleSpacing; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::Vec2ui RigContourMapGrid::numberOfElementsIJ() const +{ + return m_mapSize; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::Vec2ui RigContourMapGrid::numberOfVerticesIJ() const +{ + cvf::Vec2ui mapSize = numberOfElementsIJ(); + mapSize.x() += 1u; + mapSize.y() += 1u; + return mapSize; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +uint RigContourMapGrid::numberOfCells() const +{ + return m_mapSize.x() * m_mapSize.y(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigContourMapGrid::numberOfVertices() const +{ + cvf::Vec2ui gridSize = numberOfVerticesIJ(); + return static_cast( gridSize.x() ) * static_cast( gridSize.y() ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::Vec3d RigContourMapGrid::origin3d() const +{ + return m_expandedBoundingBox.min(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector RigContourMapGrid::generateVertices() const +{ + size_t nVertices = numberOfVertices(); + std::vector vertices( nVertices, cvf::Vec3d::ZERO ); + +#pragma omp parallel for + for ( int index = 0; index < static_cast( nVertices ); ++index ) + { + cvf::Vec2ui ij = ijFromVertexIndex( index ); + cvf::Vec2d mapPos = cellCenterPosition( ij.x(), ij.y() ); + // Shift away from sample point to vertex + mapPos.x() -= sampleSpacing() * 0.5; + mapPos.y() -= sampleSpacing() * 0.5; + + cvf::Vec3d vertexPos( mapPos, 0.0 ); + vertices[index] = vertexPos; + } + return vertices; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::Vec2ui RigContourMapGrid::ijFromVertexIndex( size_t gridIndex ) const +{ + cvf::Vec2ui gridSize = numberOfVerticesIJ(); + + uint quotientX = static_cast( gridIndex ) / gridSize.x(); + uint remainderX = static_cast( gridIndex ) % gridSize.x(); + + return cvf::Vec2ui( remainderX, quotientX ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigContourMapGrid::cellIndexFromIJ( uint i, uint j ) const +{ + CVF_ASSERT( i < m_mapSize.x() ); + CVF_ASSERT( j < m_mapSize.y() ); + + return i + j * m_mapSize.x(); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::Vec2ui RigContourMapGrid::ijFromCellIndex( size_t cellIndex ) const +{ + CVF_TIGHT_ASSERT( cellIndex < numberOfCells() ); + + uint quotientX = static_cast( cellIndex ) / m_mapSize.x(); + uint remainderX = static_cast( cellIndex ) % m_mapSize.x(); + + return cvf::Vec2ui( remainderX, quotientX ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::Vec2ui RigContourMapGrid::ijFromLocalPos( const cvf::Vec2d& localPos2d ) const +{ + uint i = localPos2d.x() / sampleSpacing(); + uint j = localPos2d.y() / sampleSpacing(); + return cvf::Vec2ui( i, j ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::Vec2d RigContourMapGrid::cellCenterPosition( uint i, uint j ) const +{ + cvf::Vec3d gridExtent = m_expandedBoundingBox.extent(); + cvf::Vec2d cellCorner = cvf::Vec2d( ( i * gridExtent.x() ) / ( m_mapSize.x() ), ( j * gridExtent.y() ) / ( m_mapSize.y() ) ); + + return cellCorner + cvf::Vec2d( sampleSpacing() * 0.5, sampleSpacing() * 0.5 ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::Vec2d RigContourMapGrid::origin2d() const +{ + return cvf::Vec2d( m_expandedBoundingBox.min().x(), m_expandedBoundingBox.min().y() ); +} + +//-------------------------------------------------------------------------------------------------- +/// Vertex positions in local coordinates (add origin2d.x() for UTM x) +//-------------------------------------------------------------------------------------------------- +std::vector RigContourMapGrid::xVertexPositions() const +{ + double gridExtent = m_expandedBoundingBox.extent().x(); + + cvf::Vec2ui gridSize = numberOfVerticesIJ(); + std::vector positions; + positions.reserve( gridSize.x() ); + for ( uint i = 0; i < gridSize.x(); ++i ) + { + positions.push_back( ( i * gridExtent ) / ( gridSize.x() - 1 ) ); + } + + return positions; +} + +//-------------------------------------------------------------------------------------------------- +/// Vertex positions in local coordinates (add origin2d.y() for UTM y) +//-------------------------------------------------------------------------------------------------- +std::vector RigContourMapGrid::yVertexPositions() const +{ + double gridExtent = m_expandedBoundingBox.extent().y(); + + cvf::Vec2ui gridSize = numberOfVerticesIJ(); + std::vector positions; + positions.reserve( gridSize.y() ); + for ( uint j = 0; j < gridSize.y(); ++j ) + { + positions.push_back( ( j * gridExtent ) / ( gridSize.y() - 1 ) ); + } + + return positions; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +double RigContourMapGrid::gridEdgeOffset() const +{ + return sampleSpacing() * 2.0; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +const cvf::BoundingBox& RigContourMapGrid::expandedBoundingBox() const +{ + return m_expandedBoundingBox; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +const cvf::BoundingBox& RigContourMapGrid::originalBoundingBox() const +{ + return m_originalBoundingBox; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +const cvf::Vec2ui& RigContourMapGrid::mapSize() const +{ + return m_mapSize; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +size_t RigContourMapGrid::vertexIndexFromIJ( uint i, uint j ) const +{ + return i + j * ( m_mapSize.x() + 1 ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::Vec2ui RigContourMapGrid::calculateMapSize( const cvf::Vec3d& gridExtent, double sampleSpacing ) +{ + uint projectionSizeX = static_cast( std::ceil( gridExtent.x() / sampleSpacing ) ); + uint projectionSizeY = static_cast( std::ceil( gridExtent.y() / sampleSpacing ) ); + + return cvf::Vec2ui( projectionSizeX, projectionSizeY ); +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +cvf::BoundingBox + RigContourMapGrid::makeMaxPointMultipleOfCellSize( const cvf::BoundingBox& boundingBox, const cvf::Vec2ui& mapSize, double sampleSpacing ) +{ + cvf::Vec3d minPoint = boundingBox.min(); + cvf::Vec3d maxPoint = boundingBox.max(); + maxPoint.x() = minPoint.x() + mapSize.x() * sampleSpacing; + maxPoint.y() = minPoint.y() + mapSize.y() * sampleSpacing; + return cvf::BoundingBox( minPoint, maxPoint ); +} diff --git a/ApplicationLibCode/ReservoirDataModel/RigContourMapGrid.h b/ApplicationLibCode/ReservoirDataModel/RigContourMapGrid.h new file mode 100644 index 0000000000..4027c5f38c --- /dev/null +++ b/ApplicationLibCode/ReservoirDataModel/RigContourMapGrid.h @@ -0,0 +1,75 @@ +///////////////////////////////////////////////////////////////////////////////// +// +// Copyright (C) 2024- Equinor ASA +// +// 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 +// for more details. +// +///////////////////////////////////////////////////////////////////////////////// + +#pragma once + +#include "cvfBoundingBox.h" +#include "cvfVector2.h" + +//================================================================================================== +/// +/// +//================================================================================================== +class RigContourMapGrid +{ +public: + RigContourMapGrid( const cvf::BoundingBox& originalBoundingBox, double sampleSpacing ); + RigContourMapGrid( const cvf::BoundingBox& originalBoundingBox, const cvf::BoundingBox& expandedBoundingBox, double sampleSpacing ); + + double sampleSpacing() const; + + cvf::Vec2ui numberOfElementsIJ() const; + cvf::Vec2ui numberOfVerticesIJ() const; + + uint numberOfCells() const; + uint numberOfValidCells() const; + size_t numberOfVertices() const; + + cvf::Vec3d origin3d() const; + + std::vector xVertexPositions() const; + std::vector yVertexPositions() const; + + std::vector generateVertices() const; + + // Cell index and position conversion + size_t cellIndexFromIJ( uint i, uint j ) const; + size_t vertexIndexFromIJ( uint i, uint j ) const; + cvf::Vec2ui ijFromVertexIndex( size_t gridIndex ) const; + cvf::Vec2ui ijFromCellIndex( size_t mapIndex ) const; + cvf::Vec2ui ijFromLocalPos( const cvf::Vec2d& localPos2d ) const; + cvf::Vec2d cellCenterPosition( uint i, uint j ) const; + cvf::Vec2d origin2d() const; + + double gridEdgeOffset() const; + + const cvf::BoundingBox& expandedBoundingBox() const; + const cvf::BoundingBox& originalBoundingBox() const; + const cvf::Vec2ui& mapSize() const; + +private: + static cvf::Vec2ui calculateMapSize( const cvf::Vec3d& extent, double sampleSpacing ); + + static cvf::BoundingBox + makeMaxPointMultipleOfCellSize( const cvf::BoundingBox& boundingBox, const cvf::Vec2ui& mapSize, double sampleSpacing ); + + double m_sampleSpacing; + cvf::Vec2ui m_mapSize; + cvf::BoundingBox m_expandedBoundingBox; + cvf::BoundingBox m_originalBoundingBox; +};