mirror of
https://github.com/OPM/ResInsight.git
synced 2025-01-01 03:37:15 -06:00
454 lines
18 KiB
C++
454 lines
18 KiB
C++
/////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Copyright (C) 2018- 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 <http://www.gnu.org/licenses/gpl.html>
|
|
// for more details.
|
|
//
|
|
/////////////////////////////////////////////////////////////////////////////////
|
|
|
|
#include "RigContourMapProjection.h"
|
|
|
|
#include "RiaStatisticsTools.h"
|
|
#include "RiaWeightedMeanCalculator.h"
|
|
|
|
#include "RigContourMapCalculator.h"
|
|
#include "RigContourMapGrid.h"
|
|
|
|
#include "cvfArray.h"
|
|
#include "cvfGeometryTools.h"
|
|
|
|
#include <algorithm>
|
|
#include <array>
|
|
#include <limits>
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
RigContourMapProjection::RigContourMapProjection( const RigContourMapGrid& contourMapGrid )
|
|
: m_contourMapGrid( contourMapGrid )
|
|
, m_currentResultTimestep( -1 )
|
|
{
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<std::vector<std::pair<size_t, double>>>
|
|
RigContourMapProjection::generateGridMapping( RigContourMapCalculator::ResultAggregationType resultAggregation,
|
|
const std::vector<double>& weights )
|
|
{
|
|
m_projected3dGridIndices = RigContourMapCalculator::generateGridMapping( *this, m_contourMapGrid, resultAggregation, weights );
|
|
return m_projected3dGridIndices;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapProjection::maxValue() const
|
|
{
|
|
return maxValue( m_aggregatedResults );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapProjection::minValue() const
|
|
{
|
|
return minValue( m_aggregatedResults );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapProjection::meanValue() const
|
|
{
|
|
return sumAllValues() / numberOfValidCells();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapProjection::sumAllValues() const
|
|
{
|
|
double sum = 0.0;
|
|
|
|
for ( size_t index = 0; index < m_aggregatedResults.size(); ++index )
|
|
{
|
|
if ( m_aggregatedResults[index] != std::numeric_limits<double>::infinity() )
|
|
{
|
|
sum += m_aggregatedResults[index];
|
|
}
|
|
}
|
|
return sum;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
cvf::Vec2ui RigContourMapProjection::numberOfElementsIJ() const
|
|
{
|
|
return m_contourMapGrid.numberOfElementsIJ();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
cvf::Vec2ui RigContourMapProjection::numberOfVerticesIJ() const
|
|
{
|
|
return m_contourMapGrid.numberOfVerticesIJ();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapProjection::valueAtVertex( unsigned int i, unsigned int j ) const
|
|
{
|
|
size_t index = m_contourMapGrid.vertexIndexFromIJ( i, j );
|
|
if ( index < numberOfVertices() )
|
|
{
|
|
return m_aggregatedVertexResults.at( index );
|
|
}
|
|
return std::numeric_limits<double>::infinity();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
unsigned int RigContourMapProjection::numberOfCells() const
|
|
{
|
|
return m_contourMapGrid.numberOfCells();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
unsigned int RigContourMapProjection::numberOfValidCells() const
|
|
{
|
|
unsigned int validCount = 0u;
|
|
for ( unsigned int i = 0; i < numberOfCells(); ++i )
|
|
{
|
|
cvf::Vec2ui ij = m_contourMapGrid.ijFromCellIndex( i );
|
|
if ( hasResultInCell( ij.x(), ij.y() ) )
|
|
{
|
|
validCount++;
|
|
}
|
|
}
|
|
return validCount;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
size_t RigContourMapProjection::numberOfVertices() const
|
|
{
|
|
cvf::Vec2ui gridSize = numberOfVerticesIJ();
|
|
return static_cast<size_t>( gridSize.x() ) * static_cast<size_t>( gridSize.y() );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool RigContourMapProjection::checkForMapIntersection( const cvf::Vec3d& domainPoint3d, cvf::Vec2d* contourMapPoint, double* valueAtPoint ) const
|
|
{
|
|
CVF_TIGHT_ASSERT( contourMapPoint );
|
|
CVF_TIGHT_ASSERT( valueAtPoint );
|
|
|
|
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<double>::infinity() )
|
|
{
|
|
*valueAtPoint = value;
|
|
*contourMapPoint = mapPos2d + gridorigin;
|
|
|
|
return true;
|
|
}
|
|
return false;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
cvf::Vec3d RigContourMapProjection::origin3d() const
|
|
{
|
|
return m_contourMapGrid.expandedBoundingBox().min();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
size_t RigContourMapProjection::gridResultIndex( size_t globalCellIdx ) const
|
|
{
|
|
return globalCellIdx;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapProjection::calculateValueInMapCell( unsigned int i,
|
|
unsigned int j,
|
|
const std::vector<double>& gridCellValues,
|
|
RigContourMapCalculator::ResultAggregationType resultAggregation ) const
|
|
{
|
|
const std::vector<std::pair<size_t, double>>& matchingCells = cellsAtIJ( i, j );
|
|
return RigContourMapCalculator::calculateValueInMapCell( *this, matchingCells, gridCellValues, resultAggregation );
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapProjection::maxValue( const std::vector<double>& aggregatedResults )
|
|
{
|
|
double maxV = RiaStatisticsTools::maximumValue( aggregatedResults );
|
|
return maxV != -std::numeric_limits<double>::max() ? maxV : -std::numeric_limits<double>::infinity();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapProjection::minValue( const std::vector<double>& aggregatedResults )
|
|
{
|
|
double minV = RiaStatisticsTools::minimumValue( aggregatedResults );
|
|
return minV != std::numeric_limits<double>::max() ? minV : std::numeric_limits<double>::infinity();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RigContourMapProjection::setCellVisibility( cvf::ref<cvf::UByteArray> cellVisibility )
|
|
{
|
|
m_cellGridIdxVisibility = cellVisibility;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
cvf::ref<cvf::UByteArray> RigContourMapProjection::getCellVisibility() const
|
|
{
|
|
return m_cellGridIdxVisibility;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RigContourMapProjection::generateVertexResults()
|
|
{
|
|
size_t nCells = numberOfCells();
|
|
if ( nCells != m_aggregatedResults.size() ) return;
|
|
|
|
size_t nVertices = numberOfVertices();
|
|
m_aggregatedVertexResults = std::vector<double>( nVertices, std::numeric_limits<double>::infinity() );
|
|
#pragma omp parallel for
|
|
for ( int index = 0; index < static_cast<int>( nVertices ); ++index )
|
|
{
|
|
cvf::Vec2ui ij = m_contourMapGrid.ijFromVertexIndex( index );
|
|
m_aggregatedVertexResults[index] = calculateValueAtVertex( ij.x(), ij.y() );
|
|
}
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapProjection::sumTriangleAreas( const std::vector<cvf::Vec4d>& triangles )
|
|
{
|
|
double sumArea = 0.0;
|
|
for ( size_t i = 0; i < triangles.size(); i += 3 )
|
|
{
|
|
cvf::Vec3d v1( triangles[i].x(), triangles[i].y(), triangles[i].z() );
|
|
cvf::Vec3d v2( triangles[i + 1].x(), triangles[i + 1].y(), triangles[i + 1].z() );
|
|
cvf::Vec3d v3( triangles[i + 2].x(), triangles[i + 2].y(), triangles[i + 2].z() );
|
|
double area = 0.5 * ( ( v3 - v1 ) ^ ( v2 - v1 ) ).length();
|
|
sumArea += area;
|
|
}
|
|
return sumArea;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapProjection::interpolateValue( const cvf::Vec2d& gridPos2d ) const
|
|
{
|
|
cvf::Vec2ui cellContainingPoint = m_contourMapGrid.ijFromLocalPos( gridPos2d );
|
|
cvf::Vec2d cellCenter = m_contourMapGrid.cellCenterPosition( cellContainingPoint.x(), cellContainingPoint.y() );
|
|
|
|
double sampleSpacing = m_contourMapGrid.sampleSpacing();
|
|
std::array<cvf::Vec3d, 4> x;
|
|
x[0] = cvf::Vec3d( cellCenter + cvf::Vec2d( -sampleSpacing * 0.5, -sampleSpacing * 0.5 ), 0.0 );
|
|
x[1] = cvf::Vec3d( cellCenter + cvf::Vec2d( sampleSpacing * 0.5, -sampleSpacing * 0.5 ), 0.0 );
|
|
x[2] = cvf::Vec3d( cellCenter + cvf::Vec2d( sampleSpacing * 0.5, sampleSpacing * 0.5 ), 0.0 );
|
|
x[3] = cvf::Vec3d( cellCenter + cvf::Vec2d( -sampleSpacing * 0.5, sampleSpacing * 0.5 ), 0.0 );
|
|
|
|
cvf::Vec4d baryCentricCoords = cvf::GeometryTools::barycentricCoords( x[0], x[1], x[2], x[3], cvf::Vec3d( gridPos2d, 0.0 ) );
|
|
|
|
std::array<cvf::Vec2ui, 4> v;
|
|
v[0] = cellContainingPoint;
|
|
v[1] = cvf::Vec2ui( cellContainingPoint.x() + 1u, cellContainingPoint.y() );
|
|
v[2] = cvf::Vec2ui( cellContainingPoint.x() + 1u, cellContainingPoint.y() + 1u );
|
|
v[3] = cvf::Vec2ui( cellContainingPoint.x(), cellContainingPoint.y() + 1u );
|
|
|
|
std::array<double, 4> vertexValues;
|
|
double validBarycentricCoordsSum = 0.0;
|
|
for ( int i = 0; i < 4; ++i )
|
|
{
|
|
double vertexValue = valueAtVertex( v[i].x(), v[i].y() );
|
|
if ( vertexValue == std::numeric_limits<double>::infinity() )
|
|
{
|
|
return std::numeric_limits<double>::infinity();
|
|
}
|
|
else
|
|
{
|
|
vertexValues[i] = vertexValue;
|
|
validBarycentricCoordsSum += baryCentricCoords[i];
|
|
}
|
|
}
|
|
|
|
if ( validBarycentricCoordsSum < 1.0e-8 )
|
|
{
|
|
return std::numeric_limits<double>::infinity();
|
|
}
|
|
|
|
// Calculate final value
|
|
double value = 0.0;
|
|
for ( int i = 0; i < 4; ++i )
|
|
{
|
|
value += baryCentricCoords[i] / validBarycentricCoordsSum * vertexValues[i];
|
|
}
|
|
|
|
return value;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapProjection::valueInCell( unsigned int i, unsigned int j ) const
|
|
{
|
|
size_t index = m_contourMapGrid.cellIndexFromIJ( i, j );
|
|
if ( index < numberOfCells() )
|
|
{
|
|
return m_aggregatedResults.at( index );
|
|
}
|
|
return std::numeric_limits<double>::infinity();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool RigContourMapProjection::hasResultInCell( unsigned int i, unsigned int j ) const
|
|
{
|
|
return !cellsAtIJ( i, j ).empty();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
double RigContourMapProjection::calculateValueAtVertex( unsigned int vi, unsigned int vj ) const
|
|
{
|
|
std::vector<unsigned int> averageIs;
|
|
std::vector<unsigned int> averageJs;
|
|
|
|
if ( vi > 0u ) averageIs.push_back( vi - 1 );
|
|
if ( vj > 0u ) averageJs.push_back( vj - 1 );
|
|
if ( vi < m_contourMapGrid.mapSize().x() ) averageIs.push_back( vi );
|
|
if ( vj < m_contourMapGrid.mapSize().y() ) averageJs.push_back( vj );
|
|
|
|
RiaWeightedMeanCalculator<double> calc;
|
|
for ( unsigned int j : averageJs )
|
|
{
|
|
for ( unsigned int i : averageIs )
|
|
{
|
|
if ( hasResultInCell( i, j ) )
|
|
{
|
|
calc.addValueAndWeight( valueInCell( i, j ), 1.0 );
|
|
}
|
|
}
|
|
}
|
|
if ( calc.validAggregatedWeight() )
|
|
{
|
|
return calc.weightedMean();
|
|
}
|
|
return std::numeric_limits<double>::infinity();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<std::pair<size_t, double>> RigContourMapProjection::cellsAtIJ( unsigned int i, unsigned int j ) const
|
|
{
|
|
size_t cellIndex = m_contourMapGrid.cellIndexFromIJ( i, j );
|
|
if ( cellIndex < m_projected3dGridIndices.size() )
|
|
{
|
|
return m_projected3dGridIndices[cellIndex];
|
|
}
|
|
return std::vector<std::pair<size_t, double>>();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// Vertex positions in local coordinates (add origin2d.x() for UTM x)
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<double> RigContourMapProjection::xVertexPositions() const
|
|
{
|
|
return m_contourMapGrid.xVertexPositions();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
/// Vertex positions in local coordinates (add origin2d.y() for UTM y)
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<double> RigContourMapProjection::yVertexPositions() const
|
|
{
|
|
return m_contourMapGrid.yVertexPositions();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
const std::vector<double>& RigContourMapProjection::aggregatedResults() const
|
|
{
|
|
return m_aggregatedResults;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
const std::vector<double>& RigContourMapProjection::aggregatedVertexResults() const
|
|
{
|
|
return m_aggregatedVertexResults;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
const std::vector<std::vector<std::pair<size_t, double>>>& RigContourMapProjection::projected3dGridIndices() const
|
|
{
|
|
return m_projected3dGridIndices;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RigContourMapProjection::clearResults()
|
|
{
|
|
m_aggregatedResults.clear();
|
|
m_aggregatedVertexResults.clear();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RigContourMapProjection::clearGridMapping()
|
|
{
|
|
m_projected3dGridIndices.clear();
|
|
}
|