mirror of
https://github.com/OPM/ResInsight.git
synced 2025-01-08 23:23:01 -06:00
504 lines
22 KiB
C++
504 lines
22 KiB
C++
/////////////////////////////////////////////////////////////////////////////////
|
|
//
|
|
// Copyright (C) Statoil 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 "RivBoxIntersectionGeometryGenerator.h"
|
|
|
|
#include "RimBoxIntersection.h"
|
|
#include "RimCase.h"
|
|
#include "RimGridView.h"
|
|
|
|
#include "cafHexGridIntersectionTools/cafHexGridIntersectionTools.h"
|
|
|
|
#include "cvfDrawableGeo.h"
|
|
#include "cvfPlane.h"
|
|
#include "cvfPrimitiveSetDirect.h"
|
|
#include "cvfStructGrid.h"
|
|
|
|
#include <array>
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
RivBoxIntersectionGeometryGenerator::RivBoxIntersectionGeometryGenerator( RimBoxIntersection* intersectionBox,
|
|
const RivIntersectionHexGridInterface* grid )
|
|
: m_intersectionBoxDefinition( intersectionBox )
|
|
, m_hexGrid( grid )
|
|
{
|
|
m_triangleVxes = new cvf::Vec3fArray;
|
|
m_cellBorderLineVxes = new cvf::Vec3fArray;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
RivBoxIntersectionGeometryGenerator::~RivBoxIntersectionGeometryGenerator()
|
|
{
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
bool RivBoxIntersectionGeometryGenerator::isAnyGeometryPresent() const
|
|
{
|
|
if ( m_triangleVxes->size() == 0 )
|
|
{
|
|
return false;
|
|
}
|
|
else
|
|
{
|
|
return true;
|
|
}
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
cvf::ref<cvf::DrawableGeo> RivBoxIntersectionGeometryGenerator::generateSurface()
|
|
{
|
|
calculateArrays();
|
|
|
|
CVF_ASSERT( m_triangleVxes.notNull() );
|
|
|
|
if ( m_triangleVxes->size() == 0 ) return nullptr;
|
|
|
|
cvf::ref<cvf::DrawableGeo> geo = new cvf::DrawableGeo;
|
|
geo->setFromTriangleVertexArray( m_triangleVxes.p() );
|
|
|
|
return geo;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
cvf::ref<cvf::DrawableGeo> RivBoxIntersectionGeometryGenerator::createMeshDrawable()
|
|
{
|
|
if ( !( m_cellBorderLineVxes.notNull() && m_cellBorderLineVxes->size() != 0 ) ) return nullptr;
|
|
|
|
cvf::ref<cvf::DrawableGeo> geo = new cvf::DrawableGeo;
|
|
geo->setVertexArray( m_cellBorderLineVxes.p() );
|
|
|
|
cvf::ref<cvf::PrimitiveSetDirect> prim = new cvf::PrimitiveSetDirect( cvf::PT_LINES );
|
|
prim->setIndexCount( m_cellBorderLineVxes->size() );
|
|
|
|
geo->addPrimitiveSet( prim.p() );
|
|
return geo;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
const std::vector<size_t>& RivBoxIntersectionGeometryGenerator::triangleToCellIndex() const
|
|
{
|
|
CVF_ASSERT( m_triangleVxes->size() );
|
|
return m_triangleToCellIdxMap;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
const std::vector<RivIntersectionVertexWeights>&
|
|
RivBoxIntersectionGeometryGenerator::triangleVxToCellCornerInterpolationWeights() const
|
|
{
|
|
CVF_ASSERT( m_triangleVxes->size() );
|
|
return m_triVxToCellCornerWeights;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
const cvf::Vec3fArray* RivBoxIntersectionGeometryGenerator::triangleVxes() const
|
|
{
|
|
CVF_ASSERT( m_triangleVxes->size() );
|
|
|
|
return m_triangleVxes.p();
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
|
|
class Box
|
|
{
|
|
public:
|
|
using FaceType = cvf::StructGridInterface::FaceType;
|
|
|
|
Box( const cvf::Mat4d& origin, const cvf::Vec3d& size )
|
|
: m_origin( origin )
|
|
, m_size( size )
|
|
{
|
|
}
|
|
|
|
std::array<cvf::Plane, 6> planes()
|
|
{
|
|
std::array<cvf::Plane, 6> boxPlanes;
|
|
|
|
boxPlanes[FaceType::POS_I].setFromPointAndNormal( m_origin.translation() + m_size, cvf::Vec3d( m_origin.col( 0 ) ) );
|
|
boxPlanes[FaceType::NEG_I].setFromPointAndNormal( m_origin.translation(), -cvf::Vec3d( m_origin.col( 0 ) ) );
|
|
boxPlanes[FaceType::POS_J].setFromPointAndNormal( m_origin.translation() + m_size, cvf::Vec3d( m_origin.col( 1 ) ) );
|
|
boxPlanes[FaceType::NEG_J].setFromPointAndNormal( m_origin.translation(), -cvf::Vec3d( m_origin.col( 1 ) ) );
|
|
boxPlanes[FaceType::POS_K].setFromPointAndNormal( m_origin.translation() + m_size, cvf::Vec3d( m_origin.col( 2 ) ) );
|
|
boxPlanes[FaceType::NEG_K].setFromPointAndNormal( m_origin.translation(), -cvf::Vec3d( m_origin.col( 2 ) ) );
|
|
|
|
return boxPlanes;
|
|
}
|
|
|
|
std::array<cvf::Vec3d, 4> faceCorners( FaceType face )
|
|
{
|
|
std::array<cvf::Vec3d, 4> corners;
|
|
|
|
cvf::Vec3d sx = cvf::Vec3d( m_origin.col( 0 ) ) * m_size[0];
|
|
cvf::Vec3d sy = cvf::Vec3d( m_origin.col( 1 ) ) * m_size[1];
|
|
cvf::Vec3d sz = cvf::Vec3d( m_origin.col( 2 ) ) * m_size[2];
|
|
|
|
switch ( face )
|
|
{
|
|
case FaceType::POS_I:
|
|
corners[0] = m_origin.translation() + sx;
|
|
corners[1] = m_origin.translation() + sx + sy;
|
|
corners[2] = m_origin.translation() + sx + sy + sz;
|
|
corners[3] = m_origin.translation() + sx + sz;
|
|
break;
|
|
case FaceType::NEG_I:
|
|
corners[0] = m_origin.translation();
|
|
corners[1] = m_origin.translation() + sz;
|
|
corners[2] = m_origin.translation() + sy + sz;
|
|
corners[3] = m_origin.translation() + sy;
|
|
break;
|
|
case FaceType::POS_J:
|
|
corners[0] = m_origin.translation() + sy;
|
|
corners[1] = m_origin.translation() + sy + sz;
|
|
corners[2] = m_origin.translation() + sy + sx + sz;
|
|
corners[3] = m_origin.translation() + sy + sx;
|
|
break;
|
|
case FaceType::NEG_J:
|
|
corners[0] = m_origin.translation();
|
|
corners[1] = m_origin.translation() + sx;
|
|
corners[2] = m_origin.translation() + sx + sz;
|
|
corners[3] = m_origin.translation() + sz;
|
|
break;
|
|
case FaceType::POS_K:
|
|
corners[0] = m_origin.translation() + sz;
|
|
corners[1] = m_origin.translation() + sz + sx;
|
|
corners[2] = m_origin.translation() + sz + sx + sy;
|
|
corners[3] = m_origin.translation() + sz + sy;
|
|
break;
|
|
case FaceType::NEG_K:
|
|
corners[0] = m_origin.translation();
|
|
corners[1] = m_origin.translation() + sy;
|
|
corners[2] = m_origin.translation() + sx + sy;
|
|
corners[3] = m_origin.translation() + sx;
|
|
break;
|
|
}
|
|
|
|
return corners;
|
|
}
|
|
|
|
// Returns the four adjacent faces in the Pos, Neg, Pos, Neg order
|
|
|
|
std::array<FaceType, 4> adjacentFaces( FaceType face )
|
|
{
|
|
std::array<FaceType, 4> clipFaces = { FaceType::NO_FACE, FaceType::NO_FACE, FaceType::NO_FACE, FaceType::NO_FACE };
|
|
FaceType oppFace = cvf::StructGridInterface::oppositeFace( face );
|
|
int clipFaceCount = 0;
|
|
for ( int faceCand = 0; faceCand < 6; ++faceCand )
|
|
{
|
|
if ( faceCand != face && faceCand != oppFace )
|
|
{
|
|
clipFaces[clipFaceCount] = (FaceType)faceCand;
|
|
clipFaceCount++;
|
|
}
|
|
}
|
|
|
|
return clipFaces;
|
|
}
|
|
|
|
private:
|
|
cvf::Mat4d m_origin;
|
|
cvf::Vec3d m_size;
|
|
};
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
RimBoxIntersection* RivBoxIntersectionGeometryGenerator::intersectionBox() const
|
|
{
|
|
return m_intersectionBoxDefinition;
|
|
}
|
|
|
|
//--------------------------------------------------------------------------------------------------
|
|
///
|
|
//--------------------------------------------------------------------------------------------------
|
|
void RivBoxIntersectionGeometryGenerator::calculateArrays()
|
|
{
|
|
if ( m_triangleVxes->size() ) return;
|
|
|
|
std::vector<cvf::Vec3f> triangleVertices;
|
|
std::vector<cvf::Vec3f> cellBorderLineVxes;
|
|
|
|
cvf::Vec3d displayOffset( 0, 0, 0 );
|
|
{
|
|
RimGridView* gridView = nullptr;
|
|
m_intersectionBoxDefinition->firstAncestorOrThisOfType( gridView );
|
|
if ( gridView && gridView->ownerCase() )
|
|
{
|
|
displayOffset = gridView->ownerCase()->displayModelOffset();
|
|
}
|
|
}
|
|
|
|
Box box( m_intersectionBoxDefinition->boxOrigin(), m_intersectionBoxDefinition->boxSize() );
|
|
std::array<cvf::Plane, 6> boxPlanes = box.planes();
|
|
|
|
RimBoxIntersection::SinglePlaneState singlePlane = m_intersectionBoxDefinition->singlePlaneState();
|
|
|
|
int startFace = 0;
|
|
int endFace = 5;
|
|
if ( singlePlane == RimBoxIntersection::PLANE_STATE_X ) startFace = endFace = Box::FaceType::POS_I;
|
|
if ( singlePlane == RimBoxIntersection::PLANE_STATE_Y ) startFace = endFace = Box::FaceType::POS_J;
|
|
if ( singlePlane == RimBoxIntersection::PLANE_STATE_Z ) startFace = endFace = Box::FaceType::POS_K;
|
|
|
|
for ( int faceIdx = startFace; faceIdx <= endFace; ++faceIdx )
|
|
{
|
|
cvf::Plane plane = boxPlanes[faceIdx];
|
|
|
|
std::array<Box::FaceType, 4> clipFaces = box.adjacentFaces( (Box::FaceType)faceIdx );
|
|
|
|
cvf::Plane p1Plane = boxPlanes[clipFaces[1]];
|
|
cvf::Plane p2Plane = boxPlanes[clipFaces[0]];
|
|
cvf::Plane p3Plane = boxPlanes[clipFaces[3]];
|
|
cvf::Plane p4Plane = boxPlanes[clipFaces[2]];
|
|
|
|
p1Plane.flip();
|
|
p2Plane.flip();
|
|
p3Plane.flip();
|
|
p4Plane.flip();
|
|
|
|
std::array<cvf::Vec3d, 4> faceCorners = box.faceCorners( (Box::FaceType)faceIdx );
|
|
|
|
cvf::BoundingBox sectionBBox;
|
|
|
|
for ( cvf::Vec3d& corner : faceCorners )
|
|
sectionBBox.add( corner );
|
|
|
|
// Similar code as IntersectionGenerator :
|
|
|
|
std::vector<size_t> columnCellCandidates;
|
|
m_hexGrid->findIntersectingCells( sectionBBox, &columnCellCandidates );
|
|
|
|
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 ( size_t cccIdx = 0; cccIdx < columnCellCandidates.size(); ++cccIdx )
|
|
{
|
|
size_t globalCellIdx = columnCellCandidates[cccIdx];
|
|
|
|
if ( !m_hexGrid->useCell( globalCellIdx ) ) continue;
|
|
|
|
hexPlaneCutTriangleVxes.clear();
|
|
m_hexGrid->cellCornerVertices( globalCellIdx, cellCorners );
|
|
m_hexGrid->cellCornerIndices( globalCellIdx, cornerIndices );
|
|
|
|
caf::HexGridIntersectionTools::planeHexIntersectionMC( plane,
|
|
cellCorners,
|
|
cornerIndices,
|
|
&hexPlaneCutTriangleVxes,
|
|
&cellFaceForEachTriangleEdge );
|
|
|
|
std::vector<caf::HexGridIntersectionTools::ClipVx> clippedTriangleVxes_once;
|
|
std::vector<int> cellFaceForEachClippedTriangleEdge_once;
|
|
caf::HexGridIntersectionTools::clipTrianglesBetweenTwoParallelPlanes( hexPlaneCutTriangleVxes,
|
|
cellFaceForEachTriangleEdge,
|
|
p1Plane,
|
|
p2Plane,
|
|
&clippedTriangleVxes_once,
|
|
&cellFaceForEachClippedTriangleEdge_once );
|
|
|
|
for ( caf::HexGridIntersectionTools::ClipVx& clvx : clippedTriangleVxes_once )
|
|
if ( !clvx.isVxIdsNative ) clvx.derivedVxLevel = 0;
|
|
|
|
std::vector<caf::HexGridIntersectionTools::ClipVx> clippedTriangleVxes;
|
|
std::vector<int> cellFaceForEachClippedTriangleEdge;
|
|
|
|
caf::HexGridIntersectionTools::clipTrianglesBetweenTwoParallelPlanes( clippedTriangleVxes_once,
|
|
cellFaceForEachClippedTriangleEdge_once,
|
|
p3Plane,
|
|
p4Plane,
|
|
&clippedTriangleVxes,
|
|
&cellFaceForEachClippedTriangleEdge );
|
|
for ( caf::HexGridIntersectionTools::ClipVx& clvx : clippedTriangleVxes )
|
|
if ( !clvx.isVxIdsNative && clvx.derivedVxLevel == -1 ) clvx.derivedVxLevel = 1;
|
|
|
|
size_t clippedTriangleCount = clippedTriangleVxes.size() / 3;
|
|
|
|
for ( uint tIdx = 0; tIdx < clippedTriangleCount; ++tIdx )
|
|
{
|
|
uint triVxIdx = tIdx * 3;
|
|
|
|
// Accumulate triangle vertices
|
|
|
|
cvf::Vec3f p0( clippedTriangleVxes[triVxIdx + 0].vx - displayOffset );
|
|
cvf::Vec3f p1( clippedTriangleVxes[triVxIdx + 1].vx - displayOffset );
|
|
cvf::Vec3f p2( clippedTriangleVxes[triVxIdx + 2].vx - displayOffset );
|
|
|
|
triangleVertices.push_back( p0 );
|
|
triangleVertices.push_back( p1 );
|
|
triangleVertices.push_back( p2 );
|
|
|
|
// Accumulate mesh lines
|
|
|
|
#define isFace( faceEnum ) ( 0 <= faceEnum && faceEnum <= 5 )
|
|
|
|
if ( isFace( cellFaceForEachClippedTriangleEdge[triVxIdx] ) )
|
|
{
|
|
cellBorderLineVxes.push_back( p0 );
|
|
cellBorderLineVxes.push_back( p1 );
|
|
}
|
|
if ( isFace( cellFaceForEachClippedTriangleEdge[triVxIdx + 1] ) )
|
|
{
|
|
cellBorderLineVxes.push_back( p1 );
|
|
cellBorderLineVxes.push_back( p2 );
|
|
}
|
|
if ( isFace( cellFaceForEachClippedTriangleEdge[triVxIdx + 2] ) )
|
|
{
|
|
cellBorderLineVxes.push_back( p2 );
|
|
cellBorderLineVxes.push_back( p0 );
|
|
}
|
|
|
|
// Mapping to cell index
|
|
|
|
m_triangleToCellIdxMap.push_back( globalCellIdx );
|
|
|
|
// Interpolation from nodes
|
|
for ( int i = 0; i < 3; ++i )
|
|
{
|
|
caf::HexGridIntersectionTools::ClipVx cvx = clippedTriangleVxes[triVxIdx + i];
|
|
if ( cvx.isVxIdsNative )
|
|
{
|
|
m_triVxToCellCornerWeights.push_back( RivIntersectionVertexWeights( cvx.clippedEdgeVx1Id,
|
|
cvx.clippedEdgeVx2Id,
|
|
cvx.normDistFromEdgeVx1 ) );
|
|
}
|
|
else
|
|
{
|
|
caf::HexGridIntersectionTools::ClipVx cvx1;
|
|
caf::HexGridIntersectionTools::ClipVx cvx2;
|
|
|
|
if ( cvx.derivedVxLevel == 0 )
|
|
{
|
|
cvx1 = hexPlaneCutTriangleVxes[cvx.clippedEdgeVx1Id];
|
|
cvx2 = hexPlaneCutTriangleVxes[cvx.clippedEdgeVx2Id];
|
|
}
|
|
else if ( cvx.derivedVxLevel == 1 )
|
|
{
|
|
cvx1 = clippedTriangleVxes_once[cvx.clippedEdgeVx1Id];
|
|
cvx2 = clippedTriangleVxes_once[cvx.clippedEdgeVx2Id];
|
|
}
|
|
else
|
|
{
|
|
CVF_ASSERT( false );
|
|
}
|
|
|
|
if ( cvx1.isVxIdsNative && cvx2.isVxIdsNative )
|
|
{
|
|
m_triVxToCellCornerWeights.push_back( RivIntersectionVertexWeights( cvx1.clippedEdgeVx1Id,
|
|
cvx1.clippedEdgeVx2Id,
|
|
cvx1.normDistFromEdgeVx1,
|
|
cvx2.clippedEdgeVx1Id,
|
|
cvx2.clippedEdgeVx2Id,
|
|
cvx2.normDistFromEdgeVx1,
|
|
cvx.normDistFromEdgeVx1 ) );
|
|
}
|
|
else
|
|
{
|
|
caf::HexGridIntersectionTools::ClipVx cvx11;
|
|
caf::HexGridIntersectionTools::ClipVx cvx12;
|
|
caf::HexGridIntersectionTools::ClipVx cvx21;
|
|
caf::HexGridIntersectionTools::ClipVx cvx22;
|
|
|
|
if ( cvx1.isVxIdsNative )
|
|
{
|
|
cvx11 = cvx1;
|
|
cvx12 = cvx1;
|
|
}
|
|
else if ( cvx1.derivedVxLevel == 0 )
|
|
{
|
|
cvx11 = hexPlaneCutTriangleVxes[cvx1.clippedEdgeVx1Id];
|
|
cvx12 = hexPlaneCutTriangleVxes[cvx1.clippedEdgeVx2Id];
|
|
}
|
|
else if ( cvx2.derivedVxLevel == 1 )
|
|
{
|
|
cvx11 = clippedTriangleVxes_once[cvx1.clippedEdgeVx1Id];
|
|
cvx12 = clippedTriangleVxes_once[cvx1.clippedEdgeVx2Id];
|
|
}
|
|
else
|
|
{
|
|
CVF_ASSERT( false );
|
|
}
|
|
|
|
if ( cvx2.isVxIdsNative )
|
|
{
|
|
cvx21 = cvx2;
|
|
cvx22 = cvx2;
|
|
}
|
|
else if ( cvx2.derivedVxLevel == 0 )
|
|
{
|
|
cvx21 = hexPlaneCutTriangleVxes[cvx2.clippedEdgeVx1Id];
|
|
cvx22 = hexPlaneCutTriangleVxes[cvx2.clippedEdgeVx2Id];
|
|
}
|
|
else if ( cvx2.derivedVxLevel == 1 )
|
|
{
|
|
cvx21 = clippedTriangleVxes_once[cvx2.clippedEdgeVx1Id];
|
|
cvx22 = clippedTriangleVxes_once[cvx2.clippedEdgeVx2Id];
|
|
}
|
|
else
|
|
{
|
|
CVF_ASSERT( false );
|
|
}
|
|
|
|
CVF_TIGHT_ASSERT( cvx11.isVxIdsNative && cvx12.isVxIdsNative && cvx21.isVxIdsNative &&
|
|
cvx22.isVxIdsNative );
|
|
|
|
m_triVxToCellCornerWeights.push_back(
|
|
RivIntersectionVertexWeights( cvx11.clippedEdgeVx1Id,
|
|
cvx11.clippedEdgeVx2Id,
|
|
cvx11.normDistFromEdgeVx1,
|
|
cvx12.clippedEdgeVx1Id,
|
|
cvx12.clippedEdgeVx2Id,
|
|
cvx2.normDistFromEdgeVx1,
|
|
cvx21.clippedEdgeVx1Id,
|
|
cvx21.clippedEdgeVx2Id,
|
|
cvx21.normDistFromEdgeVx1,
|
|
cvx22.clippedEdgeVx1Id,
|
|
cvx22.clippedEdgeVx2Id,
|
|
cvx22.normDistFromEdgeVx1,
|
|
cvx1.normDistFromEdgeVx1,
|
|
cvx2.normDistFromEdgeVx1,
|
|
cvx.normDistFromEdgeVx1 ) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
m_cellBorderLineVxes->assign( cellBorderLineVxes );
|
|
m_triangleVxes->assign( triangleVertices );
|
|
}
|