OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

Initial testing of CutAlongPloyline - WIP

Browse Source
This commit is contained in:
Jørgen Herje
2024-02-28 16:10:36 +01:00
parent 53af53024c
commit 091bd8163a
8 changed files with 522 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
ApplicationLibCode/ModelVisualization/Intersections/CMakeLists_files.cmake
View File

@@ -11,6 +11,7 @@ set(SOURCE_GROUP_HEADER_FILES
${CMAKE_CURRENT_LIST_DIR}/RivEclipseIntersectionGrid.h
${CMAKE_CURRENT_LIST_DIR}/RivFemIntersectionGrid.h
${CMAKE_CURRENT_LIST_DIR}/RivIntersectionGeometryGeneratorInterface.h
${CMAKE_CURRENT_LIST_DIR}/RivPolylineIntersectionGeometryGenerator.h
)
set(SOURCE_GROUP_SOURCE_FILES
@@ -24,6 +25,7 @@ set(SOURCE_GROUP_SOURCE_FILES
${CMAKE_CURRENT_LIST_DIR}/RivSectionFlattener.cpp
${CMAKE_CURRENT_LIST_DIR}/RivEclipseIntersectionGrid.cpp
${CMAKE_CURRENT_LIST_DIR}/RivFemIntersectionGrid.cpp
${CMAKE_CURRENT_LIST_DIR}/RivPolylineIntersectionGeometryGenerator.cpp
)
list(APPEND CODE_HEADER_FILES ${SOURCE_GROUP_HEADER_FILES})

7
ApplicationLibCode/ModelVisualization/Intersections/RivEclipseIntersectionGrid.cpp
View File

@@ -68,10 +68,15 @@ bool RivEclipseIntersectionGrid::useCell( size_t cellIndex ) const
{
const RigCell& cell = m_mainGrid->globalCellArray()[cellIndex];
if ( m_showInactiveCells )
{
return !cell.isInvalid() && ( cell.subGrid() == nullptr );
else
}
if ( m_activeCellInfo.p() != nullptr )
{
return m_activeCellInfo->isActive( cellIndex ) && ( cell.subGrid() == nullptr );
}
return true;
}
return false;
}

247
ApplicationLibCode/ModelVisualization/Intersections/RivPolylineIntersectionGeometryGenerator.cpp Normal file
View File

@@ -0,0 +1,247 @@
/////////////////////////////////////////////////////////////////////////////////
//
// 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 "RivPolylineIntersectionGeometryGenerator.h"
#include "RivIntersectionHexGridInterface.h"
#include "RivSectionFlattener.h"
#include "cafHexGridIntersectionTools/cafHexGridIntersectionTools.h"
#include "cvfPlane.h"
#pragma optimize( "", off )
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RivPolylineIntersectionGeometryGenerator::RivPolylineIntersectionGeometryGenerator( std::vector<cvf::Vec3d>& polyline,
RivIntersectionHexGridInterface* grid )
: m_polyline( polyline )
, m_hexGrid( grid )
{
m_triangleVxes = new cvf::Vec3fArray;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RivPolylineIntersectionGeometryGenerator::~RivPolylineIntersectionGeometryGenerator()
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RivPolylineIntersectionGeometryGenerator::isAnyGeometryPresent() const
{
return m_triangleVxes->size() > 0;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<size_t>& RivPolylineIntersectionGeometryGenerator::triangleToCellIndex() const
{
CVF_ASSERT( m_triangleVxes->size() > 0 );
return m_triangleToCellIdxMap;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const cvf::Vec3fArray* RivPolylineIntersectionGeometryGenerator::triangleVxes() const
{
CVF_ASSERT( m_triangleVxes->size() > 0 );
return m_triangleVxes.p();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<RivIntersectionVertexWeights>& RivPolylineIntersectionGeometryGenerator::triangleVxToCellCornerInterpolationWeights() const
{
CVF_ASSERT( m_triangleVxes->size() > 0 );
// Not implemented error
return {};
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RivPolylineIntersectionGeometryGenerator::generateIntersectionGeometry( cvf::UByteArray* visibleCells )
{
calculateArrays( visibleCells );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RivPolylineIntersectionGeometryGenerator::calculateArrays( cvf::UByteArray* visibleCells )
{
if ( m_triangleVxes->size() != 0 || m_hexGrid.isNull() ) return;
std::vector<cvf::Vec3f> calculatedTriangleVertices;
// MeshLinesAccumulator meshAcc( m_hexGrid.p() );
cvf::BoundingBox gridBBox = m_hexGrid->boundingBox();
const double topDepth = gridBBox.max().z();
const double bottomDepth = gridBBox.min().z();
std::array<cvf::Vec3d, 8> corners;
gridBBox.cornerVertices( corners.data() );
cvf::Vec3d p1_low( corners[0].x(), corners[0].y(), bottomDepth );
cvf::Vec3d p2_low( corners[1].x(), corners[1].y(), bottomDepth );
cvf::Vec3d p3_low( corners[2].x(), corners[2].y(), bottomDepth );
cvf::Plane lowPlane;
lowPlane.setFromPoints( p1_low, p2_low, p3_low );
cvf::Vec3d p1_high( p1_low.x(), p1_low.y(), topDepth );
cvf::Vec3d p2_high( p2_low.x(), p2_low.y(), topDepth );
cvf::Vec3d p3_high( p3_low.x(), p3_low.y(), topDepth );
cvf::Plane highPlane;
highPlane.setFromPoints( p1_high, p2_high, p3_high );
highPlane.flip();
const auto zAxisDirection = -cvf::Vec3d::Z_AXIS; // NOTE: Negative or positive direction?
const cvf::Vec3d maxHeightVec = zAxisDirection * gridBBox.radius();
if ( m_polyline.size() > 1 )
{
const size_t numPoints = m_polyline.size();
size_t pointIdx = 0;
while ( pointIdx < numPoints - 1 )
{
size_t nextPointIdx = RivSectionFlattener::indexToNextValidPoint( m_polyline, zAxisDirection, pointIdx );
if ( nextPointIdx == size_t( -1 ) || nextPointIdx >= m_polyline.size() )
{
break;
}
// Start and end point of polyline segment
const cvf::Vec3d p1 = m_polyline[pointIdx];
const cvf::Vec3d p2 = m_polyline[nextPointIdx];
std::vector<size_t> columnCellCandidates =
createPolylineSegmentCellCandidates( *m_hexGrid, p1, p2, maxHeightVec, topDepth, bottomDepth );
cvf::Plane plane;
plane.setFromPoints( p1, p2, p2 + maxHeightVec );
// Planes parallel to z-axis for p1 and p2, to prevent triangles outside the polyline segment
cvf::Plane p1Plane;
p1Plane.setFromPoints( p1, p1 + maxHeightVec, p1 + plane.normal() );
cvf::Plane p2Plane;
p2Plane.setFromPoints( p2, p2 + maxHeightVec, p2 - plane.normal() );
std::vector<caf::HexGridIntersectionTools::ClipVx> hexPlaneCutTriangleVxes;
hexPlaneCutTriangleVxes.reserve( 5 * 3 );
std::vector<int> cellFaceForEachTriangleEdge;
cellFaceForEachTriangleEdge.reserve( 5 * 3 );
cvf::Vec3d cellCorners[8];
size_t cornerIndices[8];
for ( const auto globalCellIdx : columnCellCandidates )
{
if ( ( visibleCells != nullptr ) && ( ( *visibleCells )[globalCellIdx] == 0 ) ) continue;
if ( !m_hexGrid->useCell( globalCellIdx ) ) continue;
hexPlaneCutTriangleVxes.clear();
m_hexGrid->cellCornerVertices( globalCellIdx, cellCorners );
m_hexGrid->cellCornerIndices( globalCellIdx, cornerIndices );
// Triangle vertices for polyline segment
caf::HexGridIntersectionTools::planeHexIntersectionMC( plane,
cellCorners,
cornerIndices,
&hexPlaneCutTriangleVxes,
&cellFaceForEachTriangleEdge );
// Clip triangles outside the polyline segment using p1 and p2 planes
std::vector<caf::HexGridIntersectionTools::ClipVx> clippedTriangleVxes;
std::vector<int> cellFaceForEachClippedTriangleEdge;
caf::HexGridIntersectionTools::clipTrianglesBetweenTwoParallelPlanes( hexPlaneCutTriangleVxes,
cellFaceForEachTriangleEdge,
p1Plane,
p2Plane,
&clippedTriangleVxes,
&cellFaceForEachClippedTriangleEdge );
for ( caf::HexGridIntersectionTools::ClipVx& clvx : clippedTriangleVxes )
{
if ( !clvx.isVxIdsNative ) clvx.derivedVxLevel = 0;
}
// Fill triangle vertices vector with clipped triangle vertices
size_t clippedTriangleCount = clippedTriangleVxes.size() / 3;
for ( size_t triangleIdx = 0; triangleIdx < clippedTriangleCount; ++triangleIdx )
{
const size_t vxIdx0 = triangleIdx * 3;
const size_t vxIdx1 = vxIdx0 + 1;
const size_t vxIdx2 = vxIdx0 + 2;
// Accumulate triangle vertices
cvf::Vec3f point0( clippedTriangleVxes[vxIdx0].vx );
cvf::Vec3f point1( clippedTriangleVxes[vxIdx1].vx );
cvf::Vec3f point2( clippedTriangleVxes[vxIdx2].vx );
calculatedTriangleVertices.emplace_back( point0 );
calculatedTriangleVertices.emplace_back( point1 );
calculatedTriangleVertices.emplace_back( point2 );
// TODO: Accumulate mesh lines?
// meshAcc.accumulateMeshLines( cellFaceForEachTriangleEdge, ...
m_triangleToCellIdxMap.push_back( globalCellIdx );
}
}
pointIdx = nextPointIdx;
}
}
m_triangleVxes->assign( calculatedTriangleVertices );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
std::vector<size_t> RivPolylineIntersectionGeometryGenerator::createPolylineSegmentCellCandidates( const RivIntersectionHexGridInterface& hexGrid,
const cvf::Vec3d& startPoint,
const cvf::Vec3d& endPoint,
const cvf::Vec3d& heightVector,
const double topDepth,
const double bottomDepth )
{
cvf::BoundingBox sectionBBox;
sectionBBox.add( startPoint );
sectionBBox.add( endPoint );
sectionBBox.add( startPoint + heightVector );
sectionBBox.add( startPoint - heightVector );
sectionBBox.add( endPoint + heightVector );
sectionBBox.add( endPoint - heightVector );
sectionBBox.cutAbove( topDepth );
sectionBBox.cutBelow( bottomDepth );
return hexGrid.findIntersectingCells( sectionBBox );
}

78
ApplicationLibCode/ModelVisualization/Intersections/RivPolylineIntersectionGeometryGenerator.h Normal file
View File

@@ -0,0 +1,78 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) Statoil ASA
// Copyright (C) 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.
//
/////////////////////////////////////////////////////////////////////////////////
#pragma once
#include "cafPdmPointer.h"
#include "RivIntersectionGeometryGeneratorInterface.h"
#include "cvfArray.h"
#include "cvfBoundingBox.h"
#include "cvfObject.h"
#include "cvfVector3.h"
#include <vector>
#include <QString>
class RigMainGrid;
class RigActiveCellInfo;
class RigResultAccessor;
class RivIntersectionHexGridInterface;
class RimSurface;
namespace cvf
{
class ScalarMapper;
class DrawableGeo;
} // namespace cvf
class RivPolylineIntersectionGeometryGenerator : public cvf::Object, public RivIntersectionGeometryGeneratorInterface
{
public:
RivPolylineIntersectionGeometryGenerator( std::vector<cvf::Vec3d>& polyline, RivIntersectionHexGridInterface* grid );
~RivPolylineIntersectionGeometryGenerator() override;
void generateIntersectionGeometry( cvf::UByteArray* visibleCells );
// GeomGen Interface
bool isAnyGeometryPresent() const override;
const std::vector<size_t>& triangleToCellIndex() const override;
const std::vector<RivIntersectionVertexWeights>& triangleVxToCellCornerInterpolationWeights() const override;
const cvf::Vec3fArray* triangleVxes() const override;
private:
void calculateArrays( cvf::UByteArray* visibleCells );
static std::vector<size_t> createPolylineSegmentCellCandidates( const RivIntersectionHexGridInterface& hexGrid,
const cvf::Vec3d& startPoint,
const cvf::Vec3d& endPoint,
const cvf::Vec3d& heightVector,
const double topDepth,
const double bottomDepth );
private:
cvf::ref<RivIntersectionHexGridInterface> m_hexGrid;
const std::vector<cvf::Vec3d> m_polyline;
// Output arrays
cvf::ref<cvf::Vec3fArray> m_triangleVxes;
std::vector<size_t> m_triangleToCellIdxMap;
};

2
ApplicationLibCode/ProjectDataModel/RimEclipseView.cpp
View File

@@ -128,6 +128,8 @@
#include <climits>
#pragma optimize( "", off )
CAF_PDM_XML_SOURCE_INIT( RimEclipseView, "ReservoirView" );
//--------------------------------------------------------------------------------------------------
///

8
GrpcInterface/GrpcProtos/GridGeometryExtraction.proto
View File

@@ -69,8 +69,14 @@ message FenceMeshSection
Vec2d endUtmXY = 7;
}
//message CutAlongPolylineResponse
//{
// repeated FenceMeshSection feceMeshSections = 1;
//GridDimensions gridDimensions = 2;
//}
message CutAlongPolylineResponse
{
repeated FenceMeshSection feceMeshSections = 1;
repeated float vertexArray = 1;
GridDimensions gridDimensions = 2;
}

82
GrpcInterface/Python/rips/WebvizPythonExamples/grid_geometry_extraction_intersection.py Normal file
View File

@@ -0,0 +1,82 @@
import numpy as np
import plotly.graph_objects as go
from rips.instance import *
from rips.generated.GridGeometryExtraction_pb2_grpc import *
from rips.generated.GridGeometryExtraction_pb2 import *
# from ..instance import *
# from ..generated.GridGeometryExtraction_pb2_grpc import *
# from ..generated.GridGeometryExtraction_pb2 import *
rips_instance = Instance.find()
grid_geometry_extraction_stub = GridGeometryExtractionStub(rips_instance.channel)
grid_file_name = None
fence_poly_line_utm_xy = [11.2631, 11.9276, 14.1083, 18.2929, 18.3523, 10.9173]
cut_along_polyline_request = GridGeometryExtraction__pb2.CutAlongPolylineRequest(
gridFilename=grid_file_name,
fencePolylineUtmXY=fence_poly_line_utm_xy,
)
cut_along_polyline_response: GridGeometryExtraction__pb2.CutAlongPolylineResponse = (
grid_geometry_extraction_stub.CutAlongPolyline(cut_along_polyline_request)
)
cut_along_polyline_response.gridDimensions
vertex_array = cut_along_polyline_response.vertexArray
num_vertex_coords = 3 # [x, y, z]
num_vertices_per_triangle = 3 # [v1, v2, v3]
num_triangles = len(vertex_array) / (num_vertex_coords * num_vertices_per_triangle)
# Create x-, y-, and z-arrays
x_array = []
y_array = []
z_array = []
for i in range(0, len(vertex_array), num_vertex_coords):
x_array.append(vertex_array[i + 0] )
y_array.append(vertex_array[i + 1] )
z_array.append(vertex_array[i + 2] )
# Create triangular mesh
i_array = []
j_array = []
k_array = []
for i in range(0, len(x_array), num_vertices_per_triangle):
# Set the indices of the vertices of the triangles
i_array.extend([i + 0])
j_array.extend([i + 1])
k_array.extend([i + 2])
fig = go.Figure(
data=[
go.Mesh3d(
x=x_array,
y=y_array,
z=z_array,
i=i_array,
j=j_array,
k=k_array,
intensity=np.linspace(-5, 5, 1000, endpoint=True),
showscale=True,
colorscale=[[0, "gold"], [0.5, "mediumturquoise"], [1.0, "magenta"]],
)
]
)
print(f"j array: {j_array}")
print(f"Number of vertices: {len(vertex_array) / 3}")
print(f"Number of traingles: {num_triangles}")
# print(f"Source cell indices array length: {len(source_cell_indices_arr)}")
# print(
# f"Origin UTM coordinates [x, y, z]: [{origin_utm.x}, {origin_utm.y}, {origin_utm.z}]"
# )
# print(
# f"Grid dimensions [I, J, K]: [{grid_dimensions.dimensions.i}, {grid_dimensions.dimensions.j}, {grid_dimensions.dimensions.k}]"
# )
print(fig.data)
fig.show()

102
GrpcInterface/RiaGrpcGridGeometryExtractionService.cpp
View File

@@ -1,6 +1,22 @@
#include "RiaGrpcGridGeometryExtractionService.h"
/////////////////////////////////////////////////////////////////////////////////
//
// 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 "VectorDefines.pb.h"
#include "RiaGrpcGridGeometryExtractionService.h"
#include "qstring.h"
@@ -9,6 +25,7 @@
#include "RiaGrpcCallbacks.h"
#include "RiaGrpcHelper.h"
#include "RifReaderSettings.h"
#include "RigActiveCellInfo.h"
#include "RigEclipseCaseData.h"
#include "RigFemPartCollection.h"
#include "RigFemPartGrid.h"
@@ -23,7 +40,9 @@
#include "RimGeoMechCase.h"
#include "RimGridView.h"
#include "RimProject.h"
#include "RivEclipseIntersectionGrid.h"
#include "RivGridPartMgr.h"
#include "RivPolylineIntersectionGeometryGenerator.h"
#include "RivReservoirPartMgr.h"
#include "RivReservoirViewPartMgr.h"
@@ -36,6 +55,8 @@
#include "qfileinfo.h"
#include "qstring.h"
#pragma optimize( "", off )
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
@@ -46,7 +67,7 @@ grpc::Status RiaGrpcGridGeometryExtractionService::GetGridSurface( grpc::ServerC
// Reset all pointers
resetInternalPointers();
// Initialize pointer
// Initialize pointers
grpc::Status status = initializeApplicationAndEclipseCaseAndEclipseViewFromAbsoluteFilePath( request->gridfilename() );
if ( status.error_code() != grpc::StatusCode::OK )
{
@@ -144,7 +165,80 @@ grpc::Status RiaGrpcGridGeometryExtractionService::CutAlongPolyline( grpc::Serve
const rips::CutAlongPolylineRequest* request,
rips::CutAlongPolylineResponse* response )
{
return grpc::Status( grpc::StatusCode::UNIMPLEMENTED, "Not implemented" );
// Reset all pointers
resetInternalPointers();
// Initialize pointers
grpc::Status status = initializeApplicationAndEclipseCaseAndEclipseViewFromAbsoluteFilePath( request->gridfilename() );
if ( status.error_code() != grpc::StatusCode::OK )
{
return status;
}
// Ensure static geometry parts created for grid part manager in view
m_eclipseView->createGridGeometryParts();
auto& fencePolyline = request->fencepolylineutmxy();
if ( fencePolyline.size() < 2 )
{
return grpc::Status( grpc::StatusCode::INVALID_ARGUMENT, "Invalid fence polyline" );
}
// Convert polyline to vector of cvf::Vec3d
std::vector<cvf::Vec3d> polyline;
const double zValue = 0.0;
for ( int i = 0; i < fencePolyline.size(); i += 2 )
{
const double xValue = fencePolyline.Get( i );
const double yValue = fencePolyline.Get( i + 1 );
cvf::Vec3d point = cvf::Vec3d( xValue, yValue, zValue );
polyline.push_back( point );
}
RigActiveCellInfo* activeCellInfo = nullptr; // No active cell info for grid
const bool showInactiveCells = true;
RivEclipseIntersectionGrid* eclipseIntersectionGrid =
new RivEclipseIntersectionGrid( m_eclipseView->mainGrid(), activeCellInfo, showInactiveCells );
auto* polylineIntersectionGenerator = new RivPolylineIntersectionGeometryGenerator( polyline, eclipseIntersectionGrid );
// Handle cell visibilities
const int firstTimeStep = 0;
cvf::UByteArray* visibleCells = new cvf::UByteArray( m_eclipseView->mainGrid()->cellCount() );
m_eclipseView->calculateCurrentTotalCellVisibility( visibleCells, firstTimeStep );
// Loop to count number of visible cells
int numVisibleCells = 0;
for ( size_t i = 0; i < visibleCells->size(); ++i )
{
if ( ( *visibleCells )[i] != 0 )
{
++numVisibleCells;
}
}
polylineIntersectionGenerator->generateIntersectionGeometry( visibleCells );
if ( !polylineIntersectionGenerator->isAnyGeometryPresent() )
{
return grpc::Status( grpc::StatusCode::INVALID_ARGUMENT, "No intersection geometry present" );
}
const auto& triangleVertices = polylineIntersectionGenerator->triangleVxes();
if ( triangleVertices->size() == 0 )
{
return grpc::Status( grpc::StatusCode::NOT_FOUND, "No triangle vertices found for polyline" );
}
// Set vertex_array and quadindicesarr response
for ( int i = 0; i < triangleVertices->size(); ++i )
{
const auto& vertex = triangleVertices->get( i );
response->add_vertexarray( vertex.x() );
response->add_vertexarray( vertex.y() );
response->add_vertexarray( vertex.z() );
}
return grpc::Status::OK;
}
//--------------------------------------------------------------------------------------------------