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

Merge pull request #5446 from OPM/feature-surface-intersection-performance

Browse Source 
Feature surface intersection performance
This commit is contained in:
Jacob Støren 2020-02-03 09:11:53 +01:00 committed by GitHub
commit 5973e8c330
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 158 additions and 113 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
ApplicationCode/ModelVisualization/Intersections/RivHexGridIntersectionTools.h
View File

@ -205,7 +205,8 @@ public:
m_weights[1] = ( (float)( normDistFromE1V1 ) );
}
explicit RivIntersectionVertexWeights( const std::array<size_t, 8> vxIds, const std::array<double, 8> explicitWeights )
explicit RivIntersectionVertexWeights( const std::array<size_t, 8>& vxIds,
const std::array<double, 8>& explicitWeights )
: m_count( 8 )
, m_vxIds( vxIds )
{

197
ApplicationCode/ModelVisualization/Surfaces/RivSurfaceIntersectionGeometryGenerator.cpp
View File

@ -144,7 +144,9 @@ void RivSurfaceIntersectionGeometryGenerator::calculateArrays()
MeshLinesAccumulator meshAcc( m_hexGrid.p() );
cvf::BoundingBox gridBBox = m_hexGrid->boundingBox();
cvf::BoundingBox gridBBox = m_hexGrid->boundingBox();
std::vector<size_t> cellDummy;
m_hexGrid->findIntersectingCells( cvf::BoundingBox( cvf::Vec3d::ZERO, cvf::Vec3d::ZERO ), &cellDummy );
m_usedSurfaceData = m_surfaceInView->surface()->surfaceData();
@ -152,103 +154,144 @@ void RivSurfaceIntersectionGeometryGenerator::calculateArrays()
const std::vector<unsigned>& nativeTriangleIndices = m_usedSurfaceData->triangleIndices();
cvf::Vec3d displayModelOffset = m_hexGrid->displayOffset();
for ( size_t ntVxIdx = 0; ntVxIdx < nativeTriangleIndices.size(); ntVxIdx += 3 )
m_triVxToCellCornerWeights.reserve( nativeTriangleIndices.size() * 24 );
outputTriangleVertices.reserve( nativeTriangleIndices.size() * 24 );
#pragma omp parallel num_threads( 6 ) // More threads have nearly no effect
{
cvf::Vec3d p0 = nativeVertices[nativeTriangleIndices[ntVxIdx + 0]];
cvf::Vec3d p1 = nativeVertices[nativeTriangleIndices[ntVxIdx + 1]];
cvf::Vec3d p2 = nativeVertices[nativeTriangleIndices[ntVxIdx + 2]];
cvf::BoundingBox triangleBBox;
triangleBBox.add( p0 );
triangleBBox.add( p1 );
triangleBBox.add( p2 );
cvf::Vec3d maxHeightVec;
std::vector<size_t> triIntersectedCellCandidates;
m_hexGrid->findIntersectingCells( triangleBBox, &triIntersectedCellCandidates );
cvf::Plane plane;
plane.setFromPoints( p0, p1, p2 );
// Loop local memory allocation.
// Must be thread private in omp parallelization
std::vector<caf::HexGridIntersectionTools::ClipVx> hexPlaneCutTriangleVxes;
hexPlaneCutTriangleVxes.reserve( 5 * 3 );
hexPlaneCutTriangleVxes.reserve( 2 * 6 * 3 );
std::vector<int> cellFaceForEachTriangleEdge;
cellFaceForEachTriangleEdge.reserve( 5 * 3 );
cellFaceForEachTriangleEdge.reserve( 2 * 6 * 3 );
std::array<cvf::Vec3d, 8> cellCorners;
std::array<size_t, 8> cornerIndices;
std::vector<cvf::Vec3d> cellCutTriangles;
cellCutTriangles.reserve( 2 * 6 * 3 );
for ( size_t ticIdx = 0; ticIdx < triIntersectedCellCandidates.size(); ++ticIdx )
std::vector<cvf::Vec3d> clippedTriangleVxes;
clippedTriangleVxes.reserve( 2 * 6 * 3 );
std::vector<int> cellFaceForEachClippedTriangleEdge;
cellFaceForEachClippedTriangleEdge.reserve( 2 * 6 * 3 );
// End loop local memory
#pragma omp for // default scheduling absolutely best
for ( int ntVxIdx = 0; ntVxIdx < nativeTriangleIndices.size(); ntVxIdx += 3 )
{
size_t globalCellIdx = triIntersectedCellCandidates[ticIdx];
cvf::Vec3d p0 = nativeVertices[nativeTriangleIndices[ntVxIdx + 0]];
cvf::Vec3d p1 = nativeVertices[nativeTriangleIndices[ntVxIdx + 1]];
cvf::Vec3d p2 = nativeVertices[nativeTriangleIndices[ntVxIdx + 2]];
if ( !m_hexGrid->useCell( globalCellIdx ) ) continue;
cvf::BoundingBox triangleBBox;
triangleBBox.add( p0 );
triangleBBox.add( p1 );
triangleBBox.add( p2 );
m_hexGrid->cellCornerVertices( globalCellIdx, &cellCorners[0] );
m_hexGrid->cellCornerIndices( globalCellIdx, &cornerIndices[0] );
cvf::Vec3d maxHeightVec;
hexPlaneCutTriangleVxes.clear();
cellFaceForEachTriangleEdge.clear();
int triangleCount = caf::HexGridIntersectionTools::planeHexIntersectionMCTet( plane,
&cellCorners[0],
&cornerIndices[0],
&hexPlaneCutTriangleVxes,
&cellFaceForEachTriangleEdge );
std::vector<size_t> triIntersectedCellCandidates;
m_hexGrid->findIntersectingCells( triangleBBox, &triIntersectedCellCandidates );
if ( triangleCount == 0 ) continue;
cvf::Plane plane;
plane.setFromPoints( p0, p1, p2 );
std::vector<cvf::Vec3d> cellCutTriangles;
for ( const auto& clipVx : hexPlaneCutTriangleVxes )
std::array<cvf::Vec3d, 8> cellCorners;
std::array<size_t, 8> cornerIndices;
for ( size_t ticIdx = 0; ticIdx < triIntersectedCellCandidates.size(); ++ticIdx )
{
cellCutTriangles.push_back( clipVx.vx );
}
size_t globalCellIdx = triIntersectedCellCandidates[ticIdx];
std::vector<cvf::Vec3d> clippedTriangleVxes;
std::vector<int> cellFaceForEachClippedTriangleEdge;
if ( !m_hexGrid->useCell( globalCellIdx ) ) continue;
caf::HexGridIntersectionTools::clipPlanarTrianglesWithInPlaneTriangle( cellCutTriangles,
cellFaceForEachTriangleEdge,
p0,
p1,
p2,
&clippedTriangleVxes,
&cellFaceForEachClippedTriangleEdge );
hexPlaneCutTriangleVxes.clear();
cellFaceForEachTriangleEdge.clear();
cellCutTriangles.clear();
size_t clippedTriangleCount = clippedTriangleVxes.size() / 3;
m_hexGrid->cellCornerVertices( globalCellIdx, &cellCorners[0] );
m_hexGrid->cellCornerIndices( globalCellIdx, &cornerIndices[0] );
for ( uint clippTrIdx = 0; clippTrIdx < clippedTriangleCount; ++clippTrIdx )
{
uint triVxIdx = clippTrIdx * 3;
int triangleCount = caf::HexGridIntersectionTools::planeHexIntersectionMCTet( plane,
&cellCorners[0],
&cornerIndices[0],
&hexPlaneCutTriangleVxes,
&cellFaceForEachTriangleEdge );
// Accumulate triangle vertices
if ( triangleCount == 0 ) continue;
cvf::Vec3d point0( clippedTriangleVxes[triVxIdx + 0] - displayModelOffset );
cvf::Vec3d point1( clippedTriangleVxes[triVxIdx + 1] - displayModelOffset );
cvf::Vec3d point2( clippedTriangleVxes[triVxIdx + 2] - displayModelOffset );
outputTriangleVertices.emplace_back( point0 );
outputTriangleVertices.emplace_back( point1 );
outputTriangleVertices.emplace_back( point2 );
// Accumulate mesh lines
meshAcc.accumulateMeshLines( cellFaceForEachClippedTriangleEdge, triVxIdx + 0, globalCellIdx, point0, point1 );
meshAcc.accumulateMeshLines( cellFaceForEachClippedTriangleEdge, triVxIdx + 1, globalCellIdx, point1, point2 );
meshAcc.accumulateMeshLines( cellFaceForEachClippedTriangleEdge, triVxIdx + 2, globalCellIdx, point2, point0 );
// Mapping to cell index
m_triangleToCellIdxMap.push_back( globalCellIdx );
// Interpolation from nodes
for ( int i = 0; i < 3; ++i )
for ( const auto& clipVx : hexPlaneCutTriangleVxes )
{
cvf::Vec3d cvx = clippedTriangleVxes[triVxIdx + i];
cellCutTriangles.push_back( clipVx.vx );
}
std::array<double, 8> cornerWeights = caf::HexInterpolator::vertexWeights( cellCorners, cvx );
m_triVxToCellCornerWeights.push_back( RivIntersectionVertexWeights( cornerIndices, cornerWeights ) );
clippedTriangleVxes.clear();
cellFaceForEachClippedTriangleEdge.clear();
caf::HexGridIntersectionTools::clipPlanarTrianglesWithInPlaneTriangle( cellCutTriangles,
cellFaceForEachTriangleEdge,
p0,
p1,
p2,
&clippedTriangleVxes,
&cellFaceForEachClippedTriangleEdge );
if ( clippedTriangleVxes.empty() ) continue;
for ( uint triVxIdx = 0; triVxIdx < clippedTriangleVxes.size(); triVxIdx += 3 )
{
// Accumulate triangle vertices
cvf::Vec3d point0( clippedTriangleVxes[triVxIdx + 0] - displayModelOffset );
cvf::Vec3d point1( clippedTriangleVxes[triVxIdx + 1] - displayModelOffset );
cvf::Vec3d point2( clippedTriangleVxes[triVxIdx + 2] - displayModelOffset );
// Interpolation from nodes
std::array<double, 8> cornerWeights0 =
caf::HexInterpolator::vertexWeights( cellCorners, clippedTriangleVxes[triVxIdx + 0] );
std::array<double, 8> cornerWeights1 =
caf::HexInterpolator::vertexWeights( cellCorners, clippedTriangleVxes[triVxIdx + 1] );
std::array<double, 8> cornerWeights2 =
caf::HexInterpolator::vertexWeights( cellCorners, clippedTriangleVxes[triVxIdx + 2] );
#pragma omp critical
{
outputTriangleVertices.emplace_back( point0 );
outputTriangleVertices.emplace_back( point1 );
outputTriangleVertices.emplace_back( point2 );
// Accumulate mesh lines
meshAcc.accumulateMeshLines( cellFaceForEachClippedTriangleEdge,
triVxIdx + 0,
globalCellIdx,
point0,
point1 );
meshAcc.accumulateMeshLines( cellFaceForEachClippedTriangleEdge,
triVxIdx + 1,
globalCellIdx,
point1,
point2 );
meshAcc.accumulateMeshLines( cellFaceForEachClippedTriangleEdge,
triVxIdx + 2,
globalCellIdx,
point2,
point0 );
// Mapping to cell index
m_triangleToCellIdxMap.push_back( globalCellIdx );
// Interpolation from nodes
m_triVxToCellCornerWeights.emplace_back(
RivIntersectionVertexWeights( cornerIndices, cornerWeights0 ) );
m_triVxToCellCornerWeights.emplace_back(
RivIntersectionVertexWeights( cornerIndices, cornerWeights1 ) );
m_triVxToCellCornerWeights.emplace_back(
RivIntersectionVertexWeights( cornerIndices, cornerWeights2 ) );
}
}
}
}

67
Fwk/AppFwk/cafHexInterpolator/cafHexInterpolator.h
View File

@ -285,6 +285,7 @@ private:
{
cvf::Vec3d Pm = interpolateInNormElm(normPoint, hexCorners);
cvf::Vec3d Pdiff = Pm - P;
if (Pdiff.lengthSquared() < 1e-3) break;
cvf::Mat3d J_inv = jacobi(hexCorners, normPoint) ;
bool inversionOk = J_inv.invert();
@ -321,44 +322,44 @@ private:
double C6_x = hexCorners[6][0]; double C6_y = hexCorners[6][1]; double C6_z = hexCorners[6][2];
double C7_x = hexCorners[7][0]; double C7_y = hexCorners[7][1]; double C7_z = hexCorners[7][2];
double dN0_di = k_1_8 * (- 1 + j + k - jk);
double dN1_di = k_1_8 * (+ 1 - j - k + jk);
double dN2_di = k_1_8 * (+ 1 + j - k - jk);
double dN3_di = k_1_8 * (- 1 - j + k + jk);
double dN4_di = k_1_8 * (- 1 + j - k + jk);
double dN5_di = k_1_8 * (+ 1 - j + k - jk);
double dN6_di = k_1_8 * (+ 1 + j + k + jk);
double dN7_di = k_1_8 * (- 1 - j - k - jk);
double dN0_di = (- 1 + j + k - jk);
double dN1_di = (+ 1 - j - k + jk);
double dN2_di = (+ 1 + j - k - jk);
double dN3_di = (- 1 - j + k + jk);
double dN4_di = (- 1 + j - k + jk);
double dN5_di = (+ 1 - j + k - jk);
double dN6_di = (+ 1 + j + k + jk);
double dN7_di = (- 1 - j - k - jk);
double dN0_dj = k_1_8 * (- 1 + i + k - ik);
double dN1_dj = k_1_8 * (- 1 - i + k + ik);
double dN2_dj = k_1_8 * (+ 1 + i - k - ik);
double dN3_dj = k_1_8 * (+ 1 - i - k + ik);
double dN4_dj = k_1_8 * (- 1 + i - k + ik);
double dN5_dj = k_1_8 * (- 1 - i - k - ik);
double dN6_dj = k_1_8 * (+ 1 + i + k + ik);
double dN7_dj = k_1_8 * (+ 1 - i + k - ik);
double dN0_dj = (- 1 + i + k - ik);
double dN1_dj = (- 1 - i + k + ik);
double dN2_dj = (+ 1 + i - k - ik);
double dN3_dj = (+ 1 - i - k + ik);
double dN4_dj = (- 1 + i - k + ik);
double dN5_dj = (- 1 - i - k - ik);
double dN6_dj = (+ 1 + i + k + ik);
double dN7_dj = (+ 1 - i + k - ik);
double dN0_dk = k_1_8 * (- 1 + i + j - ij);
double dN1_dk = k_1_8 * (- 1 - i + j + ij);
double dN2_dk = k_1_8 * (- 1 - i - j - ij);
double dN3_dk = k_1_8 * (- 1 + i - j + ij);
double dN4_dk = k_1_8 * (+ 1 - i - j + ij);
double dN5_dk = k_1_8 * (+ 1 + i - j - ij);
double dN6_dk = k_1_8 * (+ 1 + i + j + ij);
double dN7_dk = k_1_8 * (+ 1 - i + j - ij);
double dN0_dk = (- 1 + i + j - ij);
double dN1_dk = (- 1 - i + j + ij);
double dN2_dk = (- 1 - i - j - ij);
double dN3_dk = (- 1 + i - j + ij);
double dN4_dk = (+ 1 - i - j + ij);
double dN5_dk = (+ 1 + i - j - ij);
double dN6_dk = (+ 1 + i + j + ij);
double dN7_dk = (+ 1 - i + j - ij);
double dx_di = (dN0_di) * C0_x + (dN1_di) * C1_x + (dN2_di) * C2_x + (dN3_di) * C3_x + (dN4_di) * C4_x + (dN5_di) * C5_x + (dN6_di) * C6_x + (dN7_di) * C7_x;
double dx_dj = (dN0_dj) * C0_x + (dN1_dj) * C1_x + (dN2_dj) * C2_x + (dN3_dj) * C3_x + (dN4_dj) * C4_x + (dN5_dj) * C5_x + (dN6_dj) * C6_x + (dN7_dj) * C7_x;
double dx_dk = (dN0_dk) * C0_x + (dN1_dk) * C1_x + (dN2_dk) * C2_x + (dN3_dk) * C3_x + (dN4_dk) * C4_x + (dN5_dk) * C5_x + (dN6_dk) * C6_x + (dN7_dk) * C7_x;
double dx_di = k_1_8 * ( (dN0_di) * C0_x + (dN1_di) * C1_x + (dN2_di) * C2_x + (dN3_di) * C3_x + (dN4_di) * C4_x + (dN5_di) * C5_x + (dN6_di) * C6_x + (dN7_di) * C7_x );
double dx_dj = k_1_8 * ( (dN0_dj) * C0_x + (dN1_dj) * C1_x + (dN2_dj) * C2_x + (dN3_dj) * C3_x + (dN4_dj) * C4_x + (dN5_dj) * C5_x + (dN6_dj) * C6_x + (dN7_dj) * C7_x );
double dx_dk = k_1_8 * ( (dN0_dk) * C0_x + (dN1_dk) * C1_x + (dN2_dk) * C2_x + (dN3_dk) * C3_x + (dN4_dk) * C4_x + (dN5_dk) * C5_x + (dN6_dk) * C6_x + (dN7_dk) * C7_x );
double dy_di = (dN0_di) * C0_y + (dN1_di) * C1_y + (dN2_di) * C2_y + (dN3_di) * C3_y + (dN4_di) * C4_y + (dN5_di) * C5_y + (dN6_di) * C6_y + (dN7_di) * C7_y;
double dy_dj = (dN0_dj) * C0_y + (dN1_dj) * C1_y + (dN2_dj) * C2_y + (dN3_dj) * C3_y + (dN4_dj) * C4_y + (dN5_dj) * C5_y + (dN6_dj) * C6_y + (dN7_dj) * C7_y;
double dy_dk = (dN0_dk) * C0_y + (dN1_dk) * C1_y + (dN2_dk) * C2_y + (dN3_dk) * C3_y + (dN4_dk) * C4_y + (dN5_dk) * C5_y + (dN6_dk) * C6_y + (dN7_dk) * C7_y;
double dy_di = k_1_8 * ( (dN0_di) * C0_y + (dN1_di) * C1_y + (dN2_di) * C2_y + (dN3_di) * C3_y + (dN4_di) * C4_y + (dN5_di) * C5_y + (dN6_di) * C6_y + (dN7_di) * C7_y );
double dy_dj = k_1_8 * ( (dN0_dj) * C0_y + (dN1_dj) * C1_y + (dN2_dj) * C2_y + (dN3_dj) * C3_y + (dN4_dj) * C4_y + (dN5_dj) * C5_y + (dN6_dj) * C6_y + (dN7_dj) * C7_y );
double dy_dk = k_1_8 * ( (dN0_dk) * C0_y + (dN1_dk) * C1_y + (dN2_dk) * C2_y + (dN3_dk) * C3_y + (dN4_dk) * C4_y + (dN5_dk) * C5_y + (dN6_dk) * C6_y + (dN7_dk) * C7_y );
double dz_di = (dN0_di) * C0_z + (dN1_di) * C1_z + (dN2_di) * C2_z + (dN3_di) * C3_z + (dN4_di) * C4_z + (dN5_di) * C5_z + (dN6_di) * C6_z + (dN7_di) * C7_z;
double dz_dj = (dN0_dj) * C0_z + (dN1_dj) * C1_z + (dN2_dj) * C2_z + (dN3_dj) * C3_z + (dN4_dj) * C4_z + (dN5_dj) * C5_z + (dN6_dj) * C6_z + (dN7_dj) * C7_z;
double dz_dk = (dN0_dk) * C0_z + (dN1_dk) * C1_z + (dN2_dk) * C2_z + (dN3_dk) * C3_z + (dN4_dk) * C4_z + (dN5_dk) * C5_z + (dN6_dk) * C6_z + (dN7_dk) * C7_z;
double dz_di = k_1_8 * ( (dN0_di) * C0_z + (dN1_di) * C1_z + (dN2_di) * C2_z + (dN3_di) * C3_z + (dN4_di) * C4_z + (dN5_di) * C5_z + (dN6_di) * C6_z + (dN7_di) * C7_z );
double dz_dj = k_1_8 * ( (dN0_dj) * C0_z + (dN1_dj) * C1_z + (dN2_dj) * C2_z + (dN3_dj) * C3_z + (dN4_dj) * C4_z + (dN5_dj) * C5_z + (dN6_dj) * C6_z + (dN7_dj) * C7_z );
double dz_dk = k_1_8 * ( (dN0_dk) * C0_z + (dN1_dk) * C1_z + (dN2_dk) * C2_z + (dN3_dk) * C3_z + (dN4_dk) * C4_z + (dN5_dk) * C5_z + (dN6_dk) * C6_z + (dN7_dk) * C7_z );
// Do not know which ordering ends up as correct

2
Fwk/AppFwk/cafHexInterpolator/cafHexInterpolator_UnitTest/cafHexInterpolatorBasicTest.cpp
View File

@ -47,7 +47,7 @@ class HexInterpolatorTester
};
}
void testHex(const std::array<cvf::Vec3d, 8>& hexCorners, double tolerance = 1e-6)
void testHex(const std::array<cvf::Vec3d, 8>& hexCorners, double tolerance = 1e-4)
{
double result = 0.0;
result = caf::HexInterpolator::interpolateHex(hexCorners, {0.0,1.0,2.0,3.0,4.0,5.0,6.0,7.0 }, hexCorners[0]);

2
Fwk/AppFwk/cafVizExtensions/cafHexGridIntersectionTools/cafHexGridIntersectionTools.cpp
View File

@ -544,7 +544,7 @@ void HexGridIntersectionTools::clipPlanarTrianglesWithInPlaneTriangle(const std:
clipTrianglePlanes[1] = createPlaneFromEdgeAndPointInNormalDirection ( tp2, tp3, tp1 );
clipTrianglePlanes[2] = createPlaneFromEdgeAndPointInNormalDirection ( tp3, tp1, tp2 );
#define reserveSize 60
#define reserveSize 100
std::vector<cvf::Vec3d> currentInputTriangleVxes;
currentInputTriangleVxes.reserve(reserveSize);