mirror of
https://github.com/OPM/ResInsight.git
synced 2025-01-01 03:37:15 -06:00
355 lines
17 KiB
C++
355 lines
17 KiB
C++
/////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// 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 <http://www.gnu.org/licenses/gpl.html>
|
|
// for more details.
|
|
//
|
|
/////////////////////////////////////////////////////////////////////////////////
|
|
|
|
#include "RigContourMapTrianglesGenerator.h"
|
|
|
|
#include "RiaOpenMPTools.h"
|
|
|
|
#include "RigCellGeometryTools.h"
|
|
#include "RigContourMapGrid.h"
|
|
#include "RigContourMapProjection.h"
|
|
|
|
#include "cvfGeometryUtils.h"
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::vector<cvf::Vec4d>
|
|
RigContourMapTrianglesGenerator::generateTrianglesWithVertexValues( const RigContourMapGrid& contourMapGrid,
|
|
const RigContourMapProjection& contourMapProjection,
|
|
const std::vector<ContourPolygons>& contourPolygons,
|
|
const std::vector<double>& contourLevels,
|
|
const std::vector<double>& contourLevelCumulativeAreas,
|
|
bool discrete,
|
|
double sampleSpacing )
|
|
{
|
|
std::vector<cvf::Vec3d> vertices = contourMapGrid.generateVertices();
|
|
|
|
cvf::Vec2ui patchSize = contourMapGrid.numberOfVerticesIJ();
|
|
cvf::ref<cvf::UIntArray> faceList = new cvf::UIntArray;
|
|
cvf::GeometryUtils::tesselatePatchAsTriangles( patchSize.x(), patchSize.y(), 0u, true, faceList.p() );
|
|
|
|
const double cellArea = sampleSpacing * sampleSpacing;
|
|
const double areaThreshold = 1.0e-5 * 0.5 * cellArea;
|
|
|
|
std::vector<std::vector<std::vector<cvf::Vec3d>>> subtractPolygons;
|
|
if ( !contourPolygons.empty() )
|
|
{
|
|
subtractPolygons.resize( contourPolygons.size() );
|
|
for ( size_t i = 0; i < contourPolygons.size() - 1; ++i )
|
|
{
|
|
for ( size_t j = 0; j < contourPolygons[i + 1].size(); ++j )
|
|
{
|
|
subtractPolygons[i].push_back( contourPolygons[i + 1][j].vertices );
|
|
}
|
|
}
|
|
}
|
|
|
|
int numberOfThreads = RiaOpenMPTools::availableThreadCount();
|
|
|
|
std::vector<std::vector<std::vector<cvf::Vec4d>>> threadTriangles( numberOfThreads );
|
|
|
|
const std::vector<double>& aggregatedVertexResults = contourMapProjection.aggregatedVertexResults();
|
|
|
|
#pragma omp parallel
|
|
{
|
|
int myThread = RiaOpenMPTools::currentThreadIndex();
|
|
threadTriangles[myThread].resize( std::max( (size_t)1, contourPolygons.size() ) );
|
|
|
|
#pragma omp for schedule( dynamic )
|
|
for ( int64_t i = 0; i < (int64_t)faceList->size(); i += 3 )
|
|
{
|
|
std::vector<cvf::Vec3d> triangle( 3 );
|
|
std::vector<cvf::Vec4d> triangleWithValues( 3 );
|
|
bool anyValidVertex = false;
|
|
for ( size_t n = 0; n < 3; ++n )
|
|
{
|
|
uint vn = ( *faceList )[i + n];
|
|
double value = vn < aggregatedVertexResults.size() ? aggregatedVertexResults[vn] : std::numeric_limits<double>::infinity();
|
|
triangle[n] = vertices[vn];
|
|
triangleWithValues[n] = cvf::Vec4d( vertices[vn], value );
|
|
if ( value != std::numeric_limits<double>::infinity() )
|
|
{
|
|
anyValidVertex = true;
|
|
}
|
|
}
|
|
|
|
if ( !anyValidVertex )
|
|
{
|
|
continue;
|
|
}
|
|
|
|
if ( contourPolygons.empty() )
|
|
{
|
|
threadTriangles[myThread][0].insert( threadTriangles[myThread][0].end(), triangleWithValues.begin(), triangleWithValues.end() );
|
|
continue;
|
|
}
|
|
|
|
bool outsideOuterLimit = false;
|
|
for ( size_t c = 0; c < contourPolygons.size() && !outsideOuterLimit; ++c )
|
|
{
|
|
std::vector<std::vector<cvf::Vec3d>> intersectPolygons;
|
|
for ( size_t j = 0; j < contourPolygons[c].size(); ++j )
|
|
{
|
|
bool containsAtLeastOne = false;
|
|
for ( size_t t = 0; t < 3; ++t )
|
|
{
|
|
if ( contourPolygons[c][j].bbox.contains( triangle[t] ) )
|
|
{
|
|
containsAtLeastOne = true;
|
|
}
|
|
}
|
|
if ( containsAtLeastOne )
|
|
{
|
|
std::vector<std::vector<cvf::Vec3d>> clippedPolygons =
|
|
RigCellGeometryTools::intersectionWithPolygon( triangle, contourPolygons[c][j].vertices );
|
|
intersectPolygons.insert( intersectPolygons.end(), clippedPolygons.begin(), clippedPolygons.end() );
|
|
}
|
|
}
|
|
|
|
if ( intersectPolygons.empty() )
|
|
{
|
|
outsideOuterLimit = true;
|
|
continue;
|
|
}
|
|
|
|
std::vector<std::vector<cvf::Vec3d>> clippedPolygons;
|
|
|
|
if ( !subtractPolygons[c].empty() )
|
|
{
|
|
for ( const std::vector<cvf::Vec3d>& polygon : intersectPolygons )
|
|
{
|
|
std::vector<std::vector<cvf::Vec3d>> fullyClippedPolygons =
|
|
RigCellGeometryTools::subtractPolygons( polygon, subtractPolygons[c] );
|
|
clippedPolygons.insert( clippedPolygons.end(), fullyClippedPolygons.begin(), fullyClippedPolygons.end() );
|
|
}
|
|
}
|
|
else
|
|
{
|
|
clippedPolygons.swap( intersectPolygons );
|
|
}
|
|
|
|
std::vector<cvf::Vec4d> clippedTriangles;
|
|
for ( std::vector<cvf::Vec3d>& clippedPolygon : clippedPolygons )
|
|
{
|
|
std::vector<std::vector<cvf::Vec3d>> polygonTriangles;
|
|
if ( clippedPolygon.size() == 3u )
|
|
{
|
|
polygonTriangles.push_back( clippedPolygon );
|
|
}
|
|
else
|
|
{
|
|
cvf::Vec3d baryCenter = cvf::Vec3d::ZERO;
|
|
for ( size_t v = 0; v < clippedPolygon.size(); ++v )
|
|
{
|
|
cvf::Vec3d& clippedVertex = clippedPolygon[v];
|
|
baryCenter += clippedVertex;
|
|
}
|
|
baryCenter /= clippedPolygon.size();
|
|
for ( size_t v = 0; v < clippedPolygon.size(); ++v )
|
|
{
|
|
std::vector<cvf::Vec3d> clippedTriangle;
|
|
if ( v == clippedPolygon.size() - 1 )
|
|
{
|
|
clippedTriangle = { clippedPolygon[v], clippedPolygon[0], baryCenter };
|
|
}
|
|
else
|
|
{
|
|
clippedTriangle = { clippedPolygon[v], clippedPolygon[v + 1], baryCenter };
|
|
}
|
|
polygonTriangles.push_back( clippedTriangle );
|
|
}
|
|
}
|
|
for ( const std::vector<cvf::Vec3d>& polygonTriangle : polygonTriangles )
|
|
{
|
|
// Check triangle area
|
|
double area =
|
|
0.5 * ( ( polygonTriangle[1] - polygonTriangle[0] ) ^ ( polygonTriangle[2] - polygonTriangle[0] ) ).length();
|
|
if ( area < areaThreshold ) continue;
|
|
for ( const cvf::Vec3d& localVertex : polygonTriangle )
|
|
{
|
|
double value = std::numeric_limits<double>::infinity();
|
|
if ( discrete )
|
|
{
|
|
value = contourLevels[c] + 0.01 * ( contourLevels.back() - contourLevels.front() ) / contourLevels.size();
|
|
}
|
|
else
|
|
{
|
|
for ( size_t n = 0; n < 3; ++n )
|
|
{
|
|
if ( ( triangle[n] - localVertex ).length() < sampleSpacing * 0.01 &&
|
|
triangleWithValues[n].w() != std::numeric_limits<double>::infinity() )
|
|
{
|
|
value = triangleWithValues[n].w();
|
|
break;
|
|
}
|
|
}
|
|
if ( value == std::numeric_limits<double>::infinity() )
|
|
{
|
|
value = contourMapProjection.interpolateValue( cvf::Vec2d( localVertex.x(), localVertex.y() ) );
|
|
if ( value == std::numeric_limits<double>::infinity() )
|
|
{
|
|
value = contourLevels[c];
|
|
}
|
|
}
|
|
}
|
|
|
|
cvf::Vec4d globalVertex( localVertex, value );
|
|
clippedTriangles.push_back( globalVertex );
|
|
}
|
|
}
|
|
|
|
{
|
|
// Add critical section here due to a weird bug when running in a single thread
|
|
// Running multi threaded does not require this critical section, as we use a thread local data
|
|
// structure
|
|
#pragma omp critical
|
|
threadTriangles[myThread][c].insert( threadTriangles[myThread][c].end(),
|
|
clippedTriangles.begin(),
|
|
clippedTriangles.end() );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
std::vector<std::vector<cvf::Vec4d>> trianglesPerLevel( std::max( (size_t)1, contourPolygons.size() ) );
|
|
for ( size_t c = 0; c < trianglesPerLevel.size(); ++c )
|
|
{
|
|
std::vector<cvf::Vec4d> allTrianglesThisLevel;
|
|
for ( size_t i = 0; i < threadTriangles.size(); ++i )
|
|
{
|
|
allTrianglesThisLevel.insert( allTrianglesThisLevel.end(), threadTriangles[i][c].begin(), threadTriangles[i][c].end() );
|
|
}
|
|
|
|
double triangleAreasThisLevel = RigContourMapProjection::sumTriangleAreas( allTrianglesThisLevel );
|
|
if ( c >= contourLevelCumulativeAreas.size() || triangleAreasThisLevel > 1.0e-3 * contourLevelCumulativeAreas[c] )
|
|
{
|
|
trianglesPerLevel[c] = allTrianglesThisLevel;
|
|
}
|
|
}
|
|
|
|
std::vector<cvf::Vec4d> finalTriangles;
|
|
for ( size_t i = 0; i < trianglesPerLevel.size(); ++i )
|
|
{
|
|
finalTriangles.insert( finalTriangles.end(), trianglesPerLevel[i].begin(), trianglesPerLevel[i].end() );
|
|
}
|
|
|
|
return finalTriangles;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
std::pair<std::vector<RigContourMapTrianglesGenerator::ContourPolygons>, std::vector<double>>
|
|
RigContourMapTrianglesGenerator::generateContourPolygons( const RigContourMapGrid& contourMapGrid,
|
|
const RigContourMapProjection& contourMapProjection,
|
|
const std::vector<double>& initialContourLevels,
|
|
double sampleSpacing,
|
|
double sampleSpacingFactor,
|
|
bool smoothContourLines )
|
|
{
|
|
if ( contourMapProjection.minValue() != std::numeric_limits<double>::infinity() &&
|
|
contourMapProjection.maxValue() != -std::numeric_limits<double>::infinity() &&
|
|
std::fabs( contourMapProjection.maxValue() - contourMapProjection.minValue() ) > 1.0e-8 )
|
|
{
|
|
std::vector<double> contourLevels = initialContourLevels;
|
|
|
|
int nContourLevels = static_cast<int>( contourLevels.size() );
|
|
if ( nContourLevels > 2 )
|
|
{
|
|
const size_t N = contourLevels.size();
|
|
// Adjust contour levels slightly to avoid weird visual artifacts due to numerical error.
|
|
double fudgeFactor = 1.0e-3;
|
|
double fudgeAmountMin = fudgeFactor * ( contourLevels[1] - contourLevels[0] );
|
|
double fudgeAmountMax = fudgeFactor * ( contourLevels[N - 1u] - contourLevels[N - 2u] );
|
|
|
|
contourLevels.front() += fudgeAmountMin;
|
|
contourLevels.back() -= fudgeAmountMax;
|
|
|
|
double simplifyEpsilon = smoothContourLines ? 5.0e-2 * sampleSpacing : 1.0e-3 * sampleSpacing;
|
|
|
|
if ( nContourLevels >= 10 )
|
|
{
|
|
simplifyEpsilon *= 2.0;
|
|
}
|
|
if ( contourMapGrid.numberOfCells() > 100000 )
|
|
{
|
|
simplifyEpsilon *= 2.0;
|
|
}
|
|
else if ( contourMapGrid.numberOfCells() > 1000000 )
|
|
{
|
|
simplifyEpsilon *= 4.0;
|
|
}
|
|
|
|
const std::vector<double>& aggregatedVertexResults = contourMapProjection.aggregatedVertexResults();
|
|
|
|
std::vector<caf::ContourLines::ListOfLineSegments> unorderedLineSegmentsPerLevel =
|
|
caf::ContourLines::create( aggregatedVertexResults,
|
|
contourMapProjection.xVertexPositions(),
|
|
contourMapProjection.yVertexPositions(),
|
|
contourLevels );
|
|
|
|
std::vector<ContourPolygons> contourPolygons = std::vector<ContourPolygons>( unorderedLineSegmentsPerLevel.size() );
|
|
const double areaThreshold = 1.5 * ( sampleSpacing * sampleSpacing ) / ( sampleSpacingFactor * sampleSpacingFactor );
|
|
|
|
#pragma omp parallel for
|
|
for ( int i = 0; i < (int)unorderedLineSegmentsPerLevel.size(); ++i )
|
|
{
|
|
contourPolygons[i] = RigContourPolygonsTools::createContourPolygonsFromLineSegments( unorderedLineSegmentsPerLevel[i],
|
|
contourLevels[i],
|
|
areaThreshold );
|
|
|
|
if ( smoothContourLines )
|
|
{
|
|
RigContourPolygonsTools::smoothContourPolygons( contourPolygons[i], true, sampleSpacing );
|
|
}
|
|
|
|
for ( RigContourPolygonsTools::ContourPolygon& polygon : contourPolygons[i] )
|
|
{
|
|
RigCellGeometryTools::simplifyPolygon( &polygon.vertices, simplifyEpsilon );
|
|
}
|
|
}
|
|
|
|
// The clipping of contour polygons is intended to detect and fix a smoothed contour polygons
|
|
// crossing into an outer contour line. The current implementation has some side effects causing
|
|
// several contour lines to disappear. Disable this clipping for now
|
|
/*
|
|
if ( m_smoothContourLines() )
|
|
{
|
|
for ( size_t i = 1; i < contourPolygons.size(); ++i )
|
|
{
|
|
RigContourPolygonsTools::clipContourPolygons(&contourPolygons[i], &contourPolygons[i - 1] );
|
|
}
|
|
}
|
|
*/
|
|
|
|
std::vector<double> contourLevelCumulativeAreas( contourPolygons.size(), 0.0 );
|
|
for ( int64_t i = (int64_t)contourPolygons.size() - 1; i >= 0; --i )
|
|
{
|
|
double levelOuterArea = RigContourPolygonsTools::sumPolygonArea( contourPolygons[i] );
|
|
contourLevelCumulativeAreas[i] = levelOuterArea;
|
|
}
|
|
|
|
return { contourPolygons, contourLevelCumulativeAreas };
|
|
}
|
|
}
|
|
|
|
return {};
|
|
}
|