///////////////////////////////////////////////////////////////////////////////// // // 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 // for more details. // ///////////////////////////////////////////////////////////////////////////////// #include "RivBoxIntersectionGeometryGenerator.h" #include "RivIntersectionHexGridInterface.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 //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- 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 { return m_triangleVxes->size() != 0; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- cvf::ref RivBoxIntersectionGeometryGenerator::generateSurface( cvf::UByteArray* visibleCells ) { calculateArrays( visibleCells ); CVF_ASSERT( m_triangleVxes.notNull() ); if ( m_triangleVxes->size() == 0 ) return nullptr; cvf::ref geo = new cvf::DrawableGeo; geo->setFromTriangleVertexArray( m_triangleVxes.p() ); return geo; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- cvf::ref RivBoxIntersectionGeometryGenerator::createMeshDrawable() { if ( !m_cellBorderLineVxes.notNull() || m_cellBorderLineVxes->size() == 0 ) return nullptr; cvf::ref geo = new cvf::DrawableGeo; geo->setVertexArray( m_cellBorderLineVxes.p() ); cvf::ref prim = new cvf::PrimitiveSetDirect( cvf::PT_LINES ); prim->setIndexCount( m_cellBorderLineVxes->size() ); geo->addPrimitiveSet( prim.p() ); return geo; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- const std::vector& RivBoxIntersectionGeometryGenerator::triangleToCellIndex() const { CVF_ASSERT( m_triangleVxes->size() ); return m_triangleToCellIdxMap; } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- const std::vector& 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 planes() { std::array 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 faceCorners( FaceType face ) { std::array 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 adjacentFaces( FaceType face ) { std::array 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( cvf::UByteArray* visibleCells ) { if ( m_triangleVxes->size() ) return; std::vector triangleVertices; std::vector cellBorderLineVxes; cvf::Vec3d displayOffset( 0, 0, 0 ); { auto gridView = m_intersectionBoxDefinition->firstAncestorOrThisOfType(); if ( gridView && gridView->ownerCase() ) { displayOffset = gridView->ownerCase()->displayModelOffset(); } } Box box( m_intersectionBoxDefinition->boxOrigin(), m_intersectionBoxDefinition->boxSize() ); std::array 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 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 faceCorners = box.faceCorners( (Box::FaceType)faceIdx ); cvf::BoundingBox sectionBBox; for ( cvf::Vec3d& corner : faceCorners ) sectionBBox.add( corner ); // Similar code as IntersectionGenerator : std::vector columnCellCandidates = m_hexGrid->findIntersectingCells( sectionBBox ); std::vector hexPlaneCutTriangleVxes; hexPlaneCutTriangleVxes.reserve( 5 * 3 ); std::vector 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 ( ( visibleCells != nullptr ) && ( ( *visibleCells )[globalCellIdx] == 0 ) ) continue; 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 clippedTriangleVxes_once; std::vector 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 clippedTriangleVxes; std::vector 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 ); }