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

Backup - WIP

Browse Source
This commit is contained in:
Jørgen Herje
2024-03-06 15:32:23 +01:00
parent d300869217
commit 8f50da0212
11 changed files with 483 additions and 114 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
ApplicationLibCode/Commands/WellPathCommands/PointTangentManipulator/RicPointTangentManipulatorPartMgr.cpp
View File

@@ -109,7 +109,7 @@ void RicPointTangentManipulatorPartMgr::originAndTangent( cvf::Vec3d* origin, cv
//--------------------------------------------------------------------------------------------------
void RicPointTangentManipulatorPartMgr::setPolyline( const std::vector<cvf::Vec3d>& polyline )
{
m_polyline = polyline;
m_polylineUtm = polyline;
}
//--------------------------------------------------------------------------------------------------
@@ -192,10 +192,10 @@ void RicPointTangentManipulatorPartMgr::updateManipulatorFromRay( const cvf::Ray
double closestDistance = std::numeric_limits<double>::max();
cvf::Vec3d closestPoint;
for ( size_t i = 1; i < m_polyline.size(); i++ )
for ( size_t i = 1; i < m_polylineUtm.size(); i++ )
{
const auto& p1 = m_polyline[i];
const auto& p2 = m_polyline[i - 1];
const auto& p1 = m_polylineUtm[i];
const auto& p2 = m_polylineUtm[i - 1];
double normalizedIntersection;
const auto pointOnLine = cvf::GeometryTools::projectPointOnLine( p1, p2, newOrigin, &normalizedIntersection );
@@ -256,7 +256,7 @@ void RicPointTangentManipulatorPartMgr::endManipulator()
//--------------------------------------------------------------------------------------------------
void RicPointTangentManipulatorPartMgr::recreateAllGeometryAndParts()
{
if ( m_polyline.empty() )
if ( m_polylineUtm.empty() )
{
createHorizontalPlaneHandle();
createVerticalAxisHandle();
@@ -272,7 +272,7 @@ void RicPointTangentManipulatorPartMgr::recreateAllGeometryAndParts()
//--------------------------------------------------------------------------------------------------
void RicPointTangentManipulatorPartMgr::createGeometryOnly()
{
if ( m_polyline.empty() )
if ( m_polylineUtm.empty() )
{
m_handleParts[HandleType::HORIZONTAL_PLANE]->setDrawable( createHorizontalPlaneGeo().p() );
m_handleParts[HandleType::VERTICAL_AXIS]->setDrawable( createVerticalAxisGeo().p() );

2
ApplicationLibCode/Commands/WellPathCommands/PointTangentManipulator/RicPointTangentManipulatorPartMgr.h
View File

@@ -102,7 +102,7 @@ private:
double m_handleSize;
bool m_isGeometryUpdateNeeded;
std::vector<cvf::Vec3d> m_polyline;
std::vector<cvf::Vec3d> m_polylineUtm;
HandleType m_activeHandle;
cvf::Vec3d m_initialPickPoint;

31
ApplicationLibCode/ModelVisualization/Intersections/RivEnclosingPolygonGenerator.cpp
View File

@@ -25,7 +25,7 @@
///
//--------------------------------------------------------------------------------------------------
PolygonVertexWelder::PolygonVertexWelder( double weldEpsilon )
: m_epsilon( weldEpsilon )
: m_epsilonSquared( weldEpsilon * weldEpsilon )
, m_first( cvf::UNDEFINED_UINT )
{
}
@@ -52,14 +52,6 @@ std::vector<cvf::uint> PolygonVertexWelder::weldVerticesAndGetIndices( const std
//--------------------------------------------------------------------------------------------------
cvf::uint PolygonVertexWelder::weldVertexAndGetIndex( const cvf::Vec3f& vertex )
{
// Compute cell coordinates of bounding box of vertex epsilon neighborhood
int left = static_cast<int>( ( vertex.x() - m_epsilon ) );
int right = static_cast<int>( ( vertex.x() + m_epsilon ) );
int front = static_cast<int>( ( vertex.y() - m_epsilon ) );
int back = static_cast<int>( ( vertex.y() + m_epsilon ) );
int bottom = static_cast<int>( ( vertex.z() - m_epsilon ) );
int top = static_cast<int>( ( vertex.z() + m_epsilon ) );
// Call function to step through linked list of bucket, testing
// if vertex is within the epsilon of one of the vertices in the bucket
cvf::uint indexOfLocatedVertex = locateVertexInPolygon( vertex );
@@ -96,14 +88,12 @@ cvf::ref<cvf::Vec3fArray> PolygonVertexWelder::createVertexArray() const
//--------------------------------------------------------------------------------------------------
cvf::uint PolygonVertexWelder::locateVertexInPolygon( const cvf::Vec3f& vertex ) const
{
const auto epsilonSquared = m_epsilon * m_epsilon;
cvf::uint currentIndex = m_first;
while ( currentIndex != cvf::UNDEFINED_UINT )
{
// Weld point within tolerance
float distanceSquared = ( m_vertex[currentIndex] - vertex ).lengthSquared();
if ( distanceSquared < epsilonSquared )
if ( distanceSquared < m_epsilonSquared )
{
return currentIndex;
}
@@ -170,7 +160,7 @@ void RivEnclosingPolygonGenerator::constructEnclosingPolygon()
//
// Must have at least 3 edges to construct a polygon
CVF_ASSERT( m_allEdgeKeys.size() >= 3 );
CVF_ASSERT( isValidPolygon() );
// Map of edge key and number of occurrences
std::map<cvf::uint64, cvf::uint> edgeKeysAndCount;
@@ -190,13 +180,13 @@ void RivEnclosingPolygonGenerator::constructEnclosingPolygon()
}
// At least a triangle is needed to construct a polygon
CVF_ASSERT( edgeKeysAndCount.size() >= 3 ); // This occurs to often?
CVF_ASSERT( edgeKeysAndCount.size() >= 3 );
// Extract boundary edge keys from all edge keys and count
std::set<cvf::EdgeKey> boundaryEdges;
for ( const auto& [key, value] : edgeKeysAndCount )
for ( const auto& [key, count] : edgeKeysAndCount )
{
if ( value == 1 )
if ( count == 1 )
{
boundaryEdges.insert( cvf::EdgeKey::fromkeyVal( key ) );
}
@@ -220,6 +210,7 @@ void RivEnclosingPolygonGenerator::constructEnclosingPolygon()
const int end = currentEdge.index2();
if ( start == cvf::UNDEFINED_UINT || end == cvf::UNDEFINED_UINT )
{
// Throw error?
break;
}
@@ -269,6 +260,14 @@ std::vector<cvf::Vec3f> RivEnclosingPolygonGenerator::getPolygonVertices() const
return m_polygonVertices;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RivEnclosingPolygonGenerator::isValidPolygon() const
{
return m_allEdgeKeys.size() >= size_t( 3 );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

5
ApplicationLibCode/ModelVisualization/Intersections/RivEnclosingPolygonGenerator.h
View File

@@ -50,7 +50,7 @@ private:
cvf::uint addVertexToPolygon( const cvf::Vec3f& vertex );
private:
const double m_epsilon; // Tolerance for vertex welding, radius around vertex defining welding neighborhood
const double m_epsilonSquared; // Tolerance for vertex welding, radius around vertex defining welding neighborhood
cvf::uint m_first; // Start of linked list
std::vector<cvf::uint> m_next; // Links each vertex to next in linked list. Always numVertices long, will grow as vertices are added
@@ -72,8 +72,7 @@ public:
std::vector<cvf::Vec3f> getPolygonVertices() const;
bool isValidPolygon() const { return m_allEdgeKeys.size() >= size_t( 3 ); }
size_t numEdges() const { return m_allEdgeKeys.size(); }
bool isValidPolygon() const;
void addTriangleVertices( const cvf::Vec3f& p0, const cvf::Vec3f& p1, const cvf::Vec3f& p2 );
void constructEnclosingPolygon();

142
ApplicationLibCode/ModelVisualization/Intersections/RivPolylineIntersectionGeometryGenerator.cpp
View File

@@ -31,12 +31,11 @@
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RivPolylineIntersectionGeometryGenerator::RivPolylineIntersectionGeometryGenerator( std::vector<cvf::Vec3d>& polyline,
RivPolylineIntersectionGeometryGenerator::RivPolylineIntersectionGeometryGenerator( std::vector<cvf::Vec3d>& polylineUtm,
RivIntersectionHexGridInterface* grid )
: m_polyline( polyline )
: m_polylineUtm( polylineUtm )
, m_hexGrid( grid )
{
m_triangleVxes = new cvf::Vec3fArray;
m_polygonVertices = new cvf::Vec3fArray;
}
@@ -53,7 +52,6 @@ RivPolylineIntersectionGeometryGenerator::~RivPolylineIntersectionGeometryGenera
bool RivPolylineIntersectionGeometryGenerator::isAnyGeometryPresent() const
{
return m_polygonVertices->size() > 0;
// return m_triangleVxes->size() > 0;
}
//--------------------------------------------------------------------------------------------------
@@ -61,8 +59,10 @@ bool RivPolylineIntersectionGeometryGenerator::isAnyGeometryPresent() const
//--------------------------------------------------------------------------------------------------
const std::vector<size_t>& RivPolylineIntersectionGeometryGenerator::triangleToCellIndex() const
{
CVF_ASSERT( m_triangleVxes->size() > 0 );
return m_triangleToCellIdxMap;
// Not implemented - not in use
CVF_ASSERT( false );
return m_emptyTriangleToCellIdxMap;
}
//--------------------------------------------------------------------------------------------------
@@ -70,8 +70,10 @@ const std::vector<size_t>& RivPolylineIntersectionGeometryGenerator::triangleToC
//--------------------------------------------------------------------------------------------------
const cvf::Vec3fArray* RivPolylineIntersectionGeometryGenerator::triangleVxes() const
{
CVF_ASSERT( m_triangleVxes->size() > 0 );
return m_triangleVxes.p();
// Not implemented - not in use
CVF_ASSERT( false );
return nullptr;
}
//--------------------------------------------------------------------------------------------------
@@ -112,6 +114,15 @@ const std::vector<size_t>& RivPolylineIntersectionGeometryGenerator::polygonToCe
return m_polygonToCellIdxMap;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::vector<PolylineSegmentMeshData>& RivPolylineIntersectionGeometryGenerator::polylineSegmentsMeshData() const
{
CVF_ASSERT( m_polylineSegmentsMeshData.size() > 0 );
return m_polylineSegmentsMeshData;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
@@ -125,10 +136,10 @@ void RivPolylineIntersectionGeometryGenerator::generateIntersectionGeometry( cvf
//--------------------------------------------------------------------------------------------------
void RivPolylineIntersectionGeometryGenerator::calculateArrays( cvf::UByteArray* visibleCells )
{
if ( m_triangleVxes->size() != 0 || m_polygonVertices->size() != 0 || m_hexGrid.isNull() ) return;
if ( m_polygonVertices->size() != 0 || m_hexGrid.isNull() ) return;
std::vector<cvf::Vec3f> calculatedTriangleVertices;
std::vector<cvf::Vec3f> calculatedPolygonVertices;
std::vector<cvf::Vec3f> calculatedPolygonVertices;
std::vector<PolylineSegmentMeshData> polylineSegmentMeshData = {}; // Mesh data per polyline segment
cvf::BoundingBox gridBBox = m_hexGrid->boundingBox();
@@ -155,26 +166,43 @@ void RivPolylineIntersectionGeometryGenerator::calculateArrays( cvf::UByteArray*
const auto zAxisDirection = -cvf::Vec3d::Z_AXIS; // NOTE: Negative or positive direction?
const cvf::Vec3d maxHeightVec = zAxisDirection * gridBBox.radius();
if ( m_polyline.size() > 1 )
if ( m_polylineUtm.size() > 1 )
{
const size_t numPoints = m_polyline.size();
const size_t numPoints = m_polylineUtm.size();
size_t pointIdx = 0;
// Loop over polyline segments
//
// Create intersection geometry for each polyline segment and clip triangles outside the
// polyline segment, using UTM-coordinates.
//
// Afterwards convert to "local" coordinates (p1 as origin), where the local uz-coordinates
// for each vertex is calculated using Pythagoras theorem.
//
// As the plane is parallel to the z-axis, the local uz-coordinates per polyline segment
// can be calculated using Pythagoras theorem. Where a segment is defined between p1 and p2,
// which implies p1 is origin of the local coordinate system uz.
//
// For a segment a UTM coordinate (x_utm,y_utm,z_utm) converts to u = sqrt(x^2 + y^2) and z = z.
// Where x, y and z are local vertex coordinates, after subtracting the (x_utm, y_utm, z_utm)-values
// for p1 from the vertex UTM-coordinates.
while ( pointIdx < numPoints - 1 )
{
size_t nextPointIdx = RivSectionFlattener::indexToNextValidPoint( m_polyline, zAxisDirection, pointIdx );
if ( nextPointIdx == size_t( -1 ) || nextPointIdx >= m_polyline.size() )
size_t nextPointIdx = RivSectionFlattener::indexToNextValidPoint( m_polylineUtm, zAxisDirection, pointIdx );
if ( nextPointIdx == size_t( -1 ) || nextPointIdx >= m_polylineUtm.size() )
{
break;
}
// Start and end point of polyline segment
const cvf::Vec3d p1 = m_polyline[pointIdx];
const cvf::Vec3d p2 = m_polyline[nextPointIdx];
// Start and end point of polyline segment in UTM-coordinates
const cvf::Vec3d p1 = m_polylineUtm[pointIdx];
const cvf::Vec3d p2 = m_polylineUtm[nextPointIdx];
// Get cell candidates for the polyline segment (subset of global cells)
std::vector<size_t> columnCellCandidates =
createPolylineSegmentCellCandidates( *m_hexGrid, p1, p2, maxHeightVec, topDepth, bottomDepth );
// Plane for the polyline segment
cvf::Plane plane;
plane.setFromPoints( p1, p2, p2 + maxHeightVec );
@@ -192,14 +220,19 @@ void RivPolylineIntersectionGeometryGenerator::calculateArrays( cvf::UByteArray*
cvf::Vec3d cellCorners[8];
size_t cornerIndices[8];
// Mesh data for polyline segment
std::vector<cvf::Vec2f> polygonVerticesUz = {};
std::vector<cvf::uint> verticesPerPolygon = {};
const cvf::Vec3f p1_f( p1 );
// Handle triangles per cell
for ( const auto globalCellIdx : columnCellCandidates )
{
if ( ( visibleCells != nullptr ) && ( ( *visibleCells )[globalCellIdx] == 0 ) ) continue;
if ( !m_hexGrid->useCell( globalCellIdx ) ) continue;
// if ( globalCellIdx != 93996 ) continue;
// Perform intersection and clipping of triangles using UTM-coordinates
hexPlaneCutTriangleVxes.clear();
m_hexGrid->cellCornerVertices( globalCellIdx, cellCorners );
m_hexGrid->cellCornerIndices( globalCellIdx, cornerIndices );
@@ -211,18 +244,6 @@ void RivPolylineIntersectionGeometryGenerator::calculateArrays( cvf::UByteArray*
&hexPlaneCutTriangleVxes,
&cellFaceForEachTriangleEdge );
// DEBUG CODE
bool tmp = false;
if ( hexPlaneCutTriangleVxes.size() < 3 )
{
tmp = true;
}
if ( globalCellIdx == 93996 )
{
tmp = true;
}
// END DEBUG CODE
// Clip triangles outside the polyline segment using the planes for point p1 and p2
std::vector<caf::HexGridIntersectionTools::ClipVx> clippedTriangleVxes;
std::vector<int> cellFaceForEachClippedTriangleEdge;
@@ -242,7 +263,7 @@ void RivPolylineIntersectionGeometryGenerator::calculateArrays( cvf::UByteArray*
// Object to for adding triangle vertices, well vertices and generate polygon vertices
RivEnclosingPolygonGenerator enclosingPolygonGenerator;
// Fill triangle vertices vector with clipped triangle vertices
// Add clipped triangle vertices to the polygon generator using local coordinates
size_t clippedTriangleCount = clippedTriangleVxes.size() / 3;
for ( size_t triangleIdx = 0; triangleIdx < clippedTriangleCount; ++triangleIdx )
{
@@ -250,20 +271,19 @@ void RivPolylineIntersectionGeometryGenerator::calculateArrays( cvf::UByteArray*
const size_t vxIdx1 = vxIdx0 + 1;
const size_t vxIdx2 = vxIdx0 + 2;
const cvf::Vec3f point0( clippedTriangleVxes[vxIdx0].vx );
const cvf::Vec3f point1( clippedTriangleVxes[vxIdx1].vx );
const cvf::Vec3f point2( clippedTriangleVxes[vxIdx2].vx );
calculatedTriangleVertices.emplace_back( point0 );
calculatedTriangleVertices.emplace_back( point1 );
calculatedTriangleVertices.emplace_back( point2 );
// Convert points from UTM to "local" coordinates
// NOTE: Do not subtract z-values. The z-values are used for the uz-coordinates
// and not provided as additional UTM info with response
const cvf::Vec3f point0( clippedTriangleVxes[vxIdx0].vx - p1 );
const cvf::Vec3f point1( clippedTriangleVxes[vxIdx1].vx - p1 );
const cvf::Vec3f point2( clippedTriangleVxes[vxIdx2].vx - p1 );
// Add triangle to enclosing polygon line handler
enclosingPolygonGenerator.addTriangleVertices( point0, point1, point2 );
}
// Must be a triangle
if ( enclosingPolygonGenerator.numEdges() < size_t( 3 ) )
// Must be a valid polygon to continue
if ( !enclosingPolygonGenerator.isValidPolygon() )
{
continue;
}
@@ -271,24 +291,42 @@ void RivPolylineIntersectionGeometryGenerator::calculateArrays( cvf::UByteArray*
// Construct enclosing polygon after adding each triangle
enclosingPolygonGenerator.constructEnclosingPolygon();
const auto& vertices = enclosingPolygonGenerator.getPolygonVertices();
for ( const auto& vertex : enclosingPolygonGenerator.getPolygonVertices() )
{
calculatedPolygonVertices.push_back( vertex );
}
// Construct local uz-coordinates using Pythagoras theorem
for ( const auto& vertex : vertices )
{
// Convert to local uz-coordinates
const auto u = std::sqrt( vertex.x() * vertex.x() + vertex.y() * vertex.y() );
const auto z = vertex.z();
polygonVerticesUz.push_back( cvf::Vec2f( u, z ) ); // u = x, z = y
// Keep old code for debugging purposes
calculatedPolygonVertices.push_back( vertex + p1_f );
}
verticesPerPolygon.push_back( static_cast<cvf::uint>( vertices.size() ) );
// Keep old code for debugging purposes
m_verticesPerPolygon.push_back( vertices.size() );
m_polygonToCellIdxMap.push_back( globalCellIdx );
// TODO:
// - Create polygon
// - Get polygon indices
// - Convert to "local" coordinates for each polyline segment
}
// Construct polyline segment mesh data
PolylineSegmentMeshData polylineSegmentData;
polylineSegmentData.startUtmXY = cvf::Vec2d( p1.x(), p1.y() );
polylineSegmentData.endUtmXY = cvf::Vec2d( p2.x(), p2.y() );
polylineSegmentData.vertexArrayUZ.assign( polygonVerticesUz );
polylineSegmentData.verticesPerPolygon = verticesPerPolygon;
polylineSegmentData.polygonIndices = {}; // TODO: Add indices or remove from struct?
// Add polyline segment mesh data to list
m_polylineSegmentsMeshData.push_back( polylineSegmentData );
// Set next polyline segment start index
pointIdx = nextPointIdx;
}
}
m_triangleVxes->assign( calculatedTriangleVertices );
m_polygonVertices->assign( calculatedPolygonVertices );
}

24
ApplicationLibCode/ModelVisualization/Intersections/RivPolylineIntersectionGeometryGenerator.h
View File

@@ -45,10 +45,22 @@ class ScalarMapper;
class DrawableGeo;
} // namespace cvf
struct PolylineSegmentMeshData
{
// Naming things? FenceMeshSection/PolylineMeshSection?
cvf::Vec2fArray vertexArrayUZ; // Consider std::vector<cvf::uint> instead, as access methods of Vec2f is .x() and .y()
std::vector<cvf::uint> polygonIndices; // Not needed when vertices/nodes are not shared between polygons?
std::vector<cvf::uint> verticesPerPolygon;
std::vector<cvf::uint> polygonToCellIndexMap;
cvf::Vec2d startUtmXY;
cvf::Vec2d endUtmXY;
};
// TODO: Remove inheritance from RivIntersectionGeometryGeneratorInterface? As we do not use triangles
class RivPolylineIntersectionGeometryGenerator : public cvf::Object, public RivIntersectionGeometryGeneratorInterface
{
public:
RivPolylineIntersectionGeometryGenerator( std::vector<cvf::Vec3d>& polyline, RivIntersectionHexGridInterface* grid );
RivPolylineIntersectionGeometryGenerator( std::vector<cvf::Vec3d>& polylineUtm, RivIntersectionHexGridInterface* grid );
~RivPolylineIntersectionGeometryGenerator() override;
void generateIntersectionGeometry( cvf::UByteArray* visibleCells );
@@ -56,6 +68,8 @@ public:
const std::vector<size_t>& vertiesPerPolygon() const;
const std::vector<size_t>& polygonToCellIndex() const;
const std::vector<PolylineSegmentMeshData>& polylineSegmentsMeshData() const;
// GeomGen Interface
bool isAnyGeometryPresent() const override;
const std::vector<size_t>& triangleToCellIndex() const override;
@@ -73,16 +87,16 @@ private:
private:
cvf::ref<RivIntersectionHexGridInterface> m_hexGrid;
const std::vector<cvf::Vec3d> m_polyline;
const std::vector<cvf::Vec3d> m_polylineUtm;
// Output arrays
cvf::ref<cvf::Vec3fArray> m_triangleVxes;
std::vector<size_t> m_triangleToCellIdxMap;
std::vector<size_t> m_polygonToCellIdxMap;
cvf::ref<cvf::Vec3fArray> m_polygonVertices;
std::vector<size_t> m_verticesPerPolygon;
std::vector<PolylineSegmentMeshData> m_polylineSegmentsMeshData;
// Dummy vectors for GeomGen Interface
const std::vector<size_t> m_emptyTriangleToCellIdxMap = {};
const std::vector<RivIntersectionVertexWeights> m_emptyTriVxToCellCornerWeights = {};
};

19
GrpcInterface/GrpcProtos/GridGeometryExtraction.proto
View File

@@ -61,22 +61,25 @@ message FenceMeshSection
{
// U-axis defined by vector from start to end, Z is global Z
repeated float vertexArrayUZ = 1;
repeated fixed32 polyIndicesArr = 2;
repeated fixed32 polyIndicesArr = 2; // Note: Redundant if no shared nodes?
repeated fixed32 verticesPerPolygonArr = 3; // Number of vertices per polygon, numPolygons long
repeated fixed32 sourceCellIndicesArr = 4; // The originating cell index per polygon, numPolygons long
Vec2d startUtmXY = 6;
Vec2d endUtmXY = 7;
}
//message CutAlongPolylineResponse
//{
// repeated FenceMeshSection feceMeshSections = 1;
//Vec3i gridDimensions = 2;
//}
message CutAlongPolylineResponse
message PolylineTestResponse
{
// Test response returning all vertices without "local" section coordinates
repeated float polygonVertexArray = 1;
repeated fixed32 verticesPerPolygonArr = 2; // Number of vertices per polygon, numPolygons long
repeated fixed32 sourceCellIndicesArr = 3; // The originating cell index per polygon, numPolygons long
}
message CutAlongPolylineResponse
{
repeated FenceMeshSection feceMeshSections = 1;
PolylineTestResponse polylineTestResponse = 2;
Vec3i gridDimensions = 3;
}

53
GrpcInterface/Python/rips/WebvizPythonExamples/grid_geometry_extraction_cut_along_polyline.py
View File

@@ -53,9 +53,56 @@ cut_along_polyline_response: GridGeometryExtraction__pb2.CutAlongPolylineRespons
grid_geometry_extraction_stub.CutAlongPolyline(cut_along_polyline_request)
)
polygon_vertex_array_org = cut_along_polyline_response.polygonVertexArray
vertices_per_polygon = cut_along_polyline_response.verticesPerPolygonArr
source_cell_indices = cut_along_polyline_response.sourceCellIndicesArr
fence_mesh_sections = cut_along_polyline_response.feceMeshSections
print(f"Number of fence mesh sections: {len(fence_mesh_sections)}")
# for section in fence_mesh_sections:
for section in fence_mesh_sections:
polygon_vertex_array_uz = section.vertexArrayUZ
vertices_per_polygon = section.verticesPerPolygonArr
start = section.startUtmXY
end = section.endUtmXY
# Create directional vector from start to end
direction_vector = [end[0] - start[0], end[1] - start[1]]
# Decompose the polygon vertex array into x, y, z arrays
x_array = []
y_array = []
z_array = []
# 2 coordinates per vertex (u, v)
for i in range(0, len(polygon_vertex_array_uz), 2):
u = polygon_vertex_array_uz[i]
z = polygon_vertex_array_uz[i + 1]
# Calculate x, y from u and directional vector,
# where u is the length along the direction vector
x = start[0] + u * direction_vector[0]
y = start[1] + u * direction_vector[1]
x_array.append(x)
y_array.append(y)
z_array.append(z)
# ******************************************
# ******************************************
#
# TODO: CONTINUE FROM HERE
#
# ******************************************
# ******************************************
polygon_vertex_array_org = (
cut_along_polyline_response.polylineTestResponse.polygonVertexArray
)
vertices_per_polygon = (
cut_along_polyline_response.polylineTestResponse.verticesPerPolygonArr
)
source_cell_indices = (
cut_along_polyline_response.polylineTestResponse.sourceCellIndicesArr
)
x_start = polygon_vertex_array_org[0]
y_start = polygon_vertex_array_org[1]

247
GrpcInterface/Python/rips/WebvizPythonExamples/grid_geometry_extraction_cut_along_polyline_TEST_RESPONSE.py Normal file
View File

@@ -0,0 +1,247 @@
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
grid_file_name = (
"D:\\Git\\resinsight-tutorials\\model-data\\norne\\NORNE_ATW2013_RFTPLT_V2.EGRID"
)
# Test polylines
mocked_model_fence_poly_line_utm_xy = [
11.2631,
11.9276,
14.1083,
18.2929,
18.3523,
10.9173,
]
norne_case_fence_poly_line_utm_xy = [
456221,
7.32113e06,
457150,
7.32106e06,
456885,
7.32176e06,
457648,
7.3226e06,
458805,
7.32278e06,
]
norne_case_single_segment_poly_line_utm_xy = [457150, 7.32106e06, 456885, 7.32176e06]
norne_case_single_segment_poly_line_gap_utm_xy = [460877, 7.3236e06, 459279, 7.32477e06]
fence_poly_line_utm_xy = norne_case_single_segment_poly_line_utm_xy
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)
)
polygon_vertex_array_org = (
cut_along_polyline_response.polylineTestResponse.polygonVertexArray
)
vertices_per_polygon = (
cut_along_polyline_response.polylineTestResponse.verticesPerPolygonArr
)
source_cell_indices = (
cut_along_polyline_response.polylineTestResponse.sourceCellIndicesArr
)
x_start = polygon_vertex_array_org[0]
y_start = polygon_vertex_array_org[1]
z_start = polygon_vertex_array_org[2]
# Subtract x_start, y_start, z_start from all x, y, z coordinates
polygon_vertex_array = []
for i in range(0, len(polygon_vertex_array_org), 3):
polygon_vertex_array.extend(
[
polygon_vertex_array_org[i] - x_start,
polygon_vertex_array_org[i + 1] - y_start,
polygon_vertex_array_org[i + 2] - z_start,
]
)
num_vertex_coords = 3 # [x, y, z]
# Create x-, y-, and z-arrays
x_array = []
y_array = []
z_array = []
for i in range(0, len(polygon_vertex_array), num_vertex_coords):
# vertex array is provided as a single array of x, y, z coordinates
# i.e. [x1, y1, z1, x2, y2, z2, x3, y3, z3, ... , xn, yn, zn]
x_array.append(polygon_vertex_array[i + 0])
y_array.append(polygon_vertex_array[i + 1])
z_array.append(polygon_vertex_array[i + 2])
# Create triangular mesh
vertices = np.array(polygon_vertex_array).reshape(-1, 3)
# Create mesh data
x, y, z = vertices.T
i = []
j = []
k = []
# Create edges between points in triangles
triangle_edges_x = []
triangle_edges_y = []
triangle_edges_z = []
# Populate i, j, k based on vertices_per_polygon
# Create triangles from each polygon
# A quad with vertex [0,1,2,3] will be split into two triangles [0,1,2] and [0,2,3]
# A hexagon with vertex [0,1,2,3,4,5] will be split into four triangles [0,1,2], [0,2,3], [0,3,4], [0,4,5]
polygon_v0_idx = 0 # Index of vertex 0 in the polygon
for vertex_count in vertices_per_polygon:
# Must have at least one triangle
if vertex_count < 3:
polygon_v0_idx += vertex_count
continue
indices = list(range(polygon_v0_idx, polygon_v0_idx + vertex_count))
# Build triangles from polygon
num_triangles = vertex_count - 2
for triangle_index in range(0, num_triangles):
triangle_v0_idx = polygon_v0_idx
triangle_v1_idx = indices[triangle_index + 1]
triangle_v2_idx = indices[triangle_index + 2]
# Vertex indices for the triangle
i.append(triangle_v0_idx)
j.append(triangle_v1_idx)
k.append(triangle_v2_idx)
# Create edge between vertices in triangle with x,y,z coordinates, coordinates per vertex is 3
coordinate_step = 3 # step per vertex
triangle_v0_global_idx = triangle_v0_idx * coordinate_step
triangle_v1_global_idx = triangle_v1_idx * coordinate_step
triangle_v2_global_idx = triangle_v2_idx * coordinate_step
# Add x,y,z coordinates for the triangle vertices (closing triangle with 'None')
triangle_edges_x.extend(
[
polygon_vertex_array[triangle_v0_global_idx + 0],
polygon_vertex_array[triangle_v1_global_idx + 0],
polygon_vertex_array[triangle_v2_global_idx + 0],
polygon_vertex_array[triangle_v0_global_idx + 0],
None,
]
)
triangle_edges_y.extend(
[
polygon_vertex_array[triangle_v0_global_idx + 1],
polygon_vertex_array[triangle_v1_global_idx + 1],
polygon_vertex_array[triangle_v2_global_idx + 1],
polygon_vertex_array[triangle_v0_global_idx + 1],
None,
]
)
triangle_edges_z.extend(
[
polygon_vertex_array[triangle_v0_global_idx + 2],
polygon_vertex_array[triangle_v1_global_idx + 2],
polygon_vertex_array[triangle_v2_global_idx + 2],
polygon_vertex_array[triangle_v0_global_idx + 2],
None,
]
)
# Move to next polygon
polygon_v0_idx += vertex_count
# Create edges between points in polygons
polygon_edges_x = []
polygon_edges_y = []
polygon_edges_z = []
polygon_global_start_index = 0
coordinate_step = 3 # step per vertex
for vertex_count in vertices_per_polygon:
# Must have at least a triangle
if vertex_count < 3:
polygon_global_start_index += vertex_count * coordinate_step
continue
for vertex_idx in range(0, vertex_count):
vertex_global_idx = polygon_global_start_index + vertex_idx * coordinate_step
polygon_edges_x.append(polygon_vertex_array[vertex_global_idx + 0])
polygon_edges_y.append(polygon_vertex_array[vertex_global_idx + 1])
polygon_edges_z.append(polygon_vertex_array[vertex_global_idx + 2])
# Close the polygon
polygon_edges_x.append(polygon_vertex_array[polygon_global_start_index + 0])
polygon_edges_y.append(polygon_vertex_array[polygon_global_start_index + 1])
polygon_edges_z.append(polygon_vertex_array[polygon_global_start_index + 2])
polygon_edges_x.append(None)
polygon_edges_y.append(None)
polygon_edges_z.append(None)
polygon_global_start_index += vertex_count * coordinate_step
# Create mesh
mesh_3D = go.Mesh3d(
x=x, y=y, z=z, i=i, j=j, k=k, opacity=0.8, color="rgba(244,22,100,0.6)"
)
# Create edge lines for triangles
triangle_edges_3d = go.Scatter3d(
x=triangle_edges_x,
y=triangle_edges_y,
z=triangle_edges_z,
mode="lines",
name="",
line=dict(color="rgb(0,0,0)", width=1),
)
# Create outer edge lines for polygon
polygon_edges_3d = go.Scatter3d(
x=polygon_edges_x,
y=polygon_edges_y,
z=polygon_edges_z,
mode="lines",
name="",
line=dict(color="rgb(0,0,0)", width=1),
)
fig = go.Figure(
data=[
mesh_3D,
# triangle_edges_3d,
polygon_edges_3d,
]
)
# 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()

3
GrpcInterface/Python/rips/WebvizPythonExamples/grid_geometry_extraction_cut_along_polyline_TRIANGLE_VERTICES.py
View File

@@ -24,7 +24,8 @@ cut_along_polyline_response: GridGeometryExtraction__pb2.CutAlongPolylineRespons
grid_geometry_extraction_stub.CutAlongPolyline(cut_along_polyline_request)
)
vertex_array = cut_along_polyline_response.triangleVertexArray
vertex_array = cut_along_polyline_response.polylineTestResponse.triangleVertexArray
cut_along_polyline_response.polylineTestResponse
num_vertex_coords = 3 # [x, y, z]
num_vertices_per_triangle = 3 # [v1, v2, v3]

59
GrpcInterface/RiaGrpcGridGeometryExtractionService.cpp
View File

@@ -220,21 +220,40 @@ grpc::Status RiaGrpcGridGeometryExtractionService::CutAlongPolyline( grpc::Serve
return grpc::Status( grpc::StatusCode::INVALID_ARGUMENT, "No intersection geometry present" );
}
// Get results
// const auto& triangleVertices = polylineIntersectionGenerator->triangleVxes();
// if ( triangleVertices->size() == 0 )
//{
// return grpc::Status( grpc::StatusCode::NOT_FOUND, "No triangle vertices found for polyline" );
//}
// Add fence mesh sections
const auto& polylineSegmentsMeshData = polylineIntersectionGenerator->polylineSegmentsMeshData();
for ( const auto& segment : polylineSegmentsMeshData )
{
auto* fenceMeshSection = response->add_fecemeshsections();
//// Set vertex_array and quadindicesarr response
// for ( int i = 0; i < triangleVertices->size(); ++i )
//{
// const auto& vertex = triangleVertices->get( i );
// response->add_trianglevertexarray( vertex.x() );
// response->add_trianglevertexarray( vertex.y() );
// response->add_trianglevertexarray( vertex.z() );
// }
// Set start UTM (x,y)
rips::Vec2d* startUtmXY = new rips::Vec2d;
startUtmXY->set_x( segment.startUtmXY.x() );
startUtmXY->set_y( segment.startUtmXY.y() );
fenceMeshSection->set_allocated_startutmxy( startUtmXY );
// Set end UTM (x,y)
rips::Vec2d* endUtmXY = new rips::Vec2d;
endUtmXY->set_x( segment.endUtmXY.x() );
endUtmXY->set_y( segment.endUtmXY.y() );
fenceMeshSection->set_allocated_endutmxy( endUtmXY );
// Fill the vertext array
for ( const auto& vertex : segment.vertexArrayUZ )
{
fenceMeshSection->add_vertexarrayuz( vertex.x() );
fenceMeshSection->add_vertexarrayuz( vertex.y() );
}
// Fill the source cell indices array
for ( const auto& sourceCellIndex : segment.polygonToCellIndexMap )
{
fenceMeshSection->add_sourcecellindicesarr( static_cast<google::protobuf::uint32>( sourceCellIndex ) );
}
}
// Add test response
rips::PolylineTestResponse* polylineTestResponse = new rips::PolylineTestResponse;
// Polygon vertices
const auto& polygonVertices = polylineIntersectionGenerator->polygonVxes();
@@ -245,25 +264,27 @@ grpc::Status RiaGrpcGridGeometryExtractionService::CutAlongPolyline( grpc::Serve
for ( int i = 0; i < polygonVertices->size(); ++i )
{
const auto& vertex = polygonVertices->get( i );
response->add_polygonvertexarray( vertex.x() );
response->add_polygonvertexarray( vertex.y() );
response->add_polygonvertexarray( vertex.z() );
polylineTestResponse->add_polygonvertexarray( vertex.x() );
polylineTestResponse->add_polygonvertexarray( vertex.y() );
polylineTestResponse->add_polygonvertexarray( vertex.z() );
}
// Vertices per polygon
const auto& verticesPerPolygon = polylineIntersectionGenerator->vertiesPerPolygon();
for ( const auto& elm : verticesPerPolygon )
{
response->add_verticesperpolygonarr( static_cast<google::protobuf::uint32>( elm ) );
polylineTestResponse->add_verticesperpolygonarr( static_cast<google::protobuf::uint32>( elm ) );
}
// Polygon to cell indices
const auto& polygonCellIndices = polylineIntersectionGenerator->polygonToCellIndex();
for ( const auto& elm : polygonCellIndices )
{
response->add_sourcecellindicesarr( static_cast<google::protobuf::uint32>( elm ) );
polylineTestResponse->add_sourcecellindicesarr( static_cast<google::protobuf::uint32>( elm ) );
}
response->set_allocated_polylinetestresponse( polylineTestResponse );
return grpc::Status::OK;
}