From 3c405528bbf4a0d2edb77ac079306f8cb78895bf Mon Sep 17 00:00:00 2001 From: jonjenssen Date: Tue, 27 Feb 2024 01:22:25 +0100 Subject: [PATCH] Work in progress --- .../RifFaultReactivationModelExporter.cpp | 8 +-- .../Faults/RimFaultReactivationEnums.h | 3 +- .../RigFaultReactivationModel.cpp | 9 ++- .../RigFaultReactivationModelGenerator.cpp | 6 +- .../ReservoirDataModel/RigGriddedPart3d.cpp | 60 ++++++++++++++++--- .../ReservoirDataModel/RigGriddedPart3d.h | 6 +- Fwk/AppFwk/cafVizExtensions/cafLine.h | 2 + Fwk/AppFwk/cafVizExtensions/cafLine.inl | 24 ++++++++ 8 files changed, 100 insertions(+), 18 deletions(-) diff --git a/ApplicationLibCode/FileInterface/RifFaultReactivationModelExporter.cpp b/ApplicationLibCode/FileInterface/RifFaultReactivationModelExporter.cpp index 3f97281136..1f3c59557a 100644 --- a/ApplicationLibCode/FileInterface/RifFaultReactivationModelExporter.cpp +++ b/ApplicationLibCode/FileInterface/RifFaultReactivationModelExporter.cpp @@ -204,13 +204,13 @@ std::pair RifFaultReactivationModelExporter::printParts( for ( auto [boundary, boundaryName] : boundaries ) { // Create boundary condition sets for each side of the parts (except top). - auto boundaryNodes = grid->boundaryNodes(); - auto boundaryElements = grid->boundaryElements(); + const auto& boundaryNodes = grid->boundaryNodes(); + const auto& boundaryElements = grid->boundaryElements(); - const std::vector& nodes = boundaryNodes[boundary]; + const std::vector& nodes = boundaryNodes.at( boundary ); RifInpExportTools::printNodeSet( stream, boundaryName, false, nodes ); - const std::vector& elements = boundaryElements[boundary]; + const std::vector& elements = boundaryElements.at( boundary ); RifInpExportTools::printElementSet( stream, boundaryName, false, elements ); } diff --git a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationEnums.h b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationEnums.h index cb2d84e242..6bc81eef46 100644 --- a/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationEnums.h +++ b/ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationEnums.h @@ -39,7 +39,8 @@ enum class Boundary { FarSide, Bottom, - Fault + Fault, + Reservoir }; enum class ElementSets diff --git a/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.cpp b/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.cpp index 286e91b986..3d1e4d7e76 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.cpp @@ -25,6 +25,8 @@ #include "RimEclipseCase.h" +#include "cafLine.h" + #include //-------------------------------------------------------------------------------------------------- @@ -259,8 +261,13 @@ void RigFaultReactivationModel::postProcessElementSets( const RimEclipseCase* eC auto cellInfo = eCase->eclipseCaseData()->activeCellInfo( RiaDefines::PorosityModelType::MATRIX_MODEL ); + auto [top, bottom] = faultTopBottom(); + caf::Line line( top, bottom ); + for ( auto part : allGridParts() ) { - m_3dparts[part]->postProcessElementSets( eCase->mainGrid(), cellInfo ); + auto gridPart = m_3dparts[part]; + gridPart->postProcessElementSets( eCase->mainGrid(), cellInfo ); + // gridPart->postProcessBoundaryNodes( line ); } } diff --git a/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp b/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp index aac975721e..26a3faabec 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp @@ -582,7 +582,8 @@ void RigFaultReactivationModelGenerator::generateGeometry( size_t sta m_modelThickness, m_topReservoirFront.z(), m_modelNormal, - m_faultZoneCells ); + m_faultZoneCells, + faultTopBottomPoints() ); backPart->generateGeometry( m_backPoints, backReservoirLayers, m_maxCellHeight, @@ -591,7 +592,8 @@ void RigFaultReactivationModelGenerator::generateGeometry( size_t sta m_modelThickness, m_topReservoirBack.z(), -1.0 * m_modelNormal, - m_faultZoneCells ); + m_faultZoneCells, + faultTopBottomPoints() ); frontPart->generateLocalNodes( m_localCoordTransform ); backPart->generateLocalNodes( m_localCoordTransform ); diff --git a/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp b/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp index f98e3c19f8..6d5af77743 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp @@ -27,6 +27,8 @@ #include "cvfBoundingBox.h" #include "cvfTextureImage.h" +#include "cafLine.h" + #include #include @@ -207,7 +209,8 @@ void RigGriddedPart3d::generateGeometry( const std::array& input double modelThickness, double topHeight, cvf::Vec3d thicknessDirection, - int nFaultZoneCells ) + int nFaultZoneCells, + std::pair topBottomFaultPoints ) { reset(); @@ -224,9 +227,10 @@ void RigGriddedPart3d::generateGeometry( const std::array& input layersPerRegion[Regions::LowerUnderburden].pop_back(); // to avoid overlap with bottom of next region layersPerRegion[Regions::Reservoir].pop_back(); // to avoid overlap with bottom of next region - m_boundaryNodes[Boundary::Bottom] = {}; - m_boundaryNodes[Boundary::FarSide] = {}; - m_boundaryNodes[Boundary::Fault] = {}; + m_boundaryNodes[Boundary::Bottom] = {}; + m_boundaryNodes[Boundary::FarSide] = {}; + m_boundaryNodes[Boundary::Fault] = {}; + m_boundaryNodes[Boundary::Reservoir] = {}; size_t nVertCells = 0; const size_t nHorzCells = horizontalPartition.size() - 1; @@ -236,9 +240,16 @@ void RigGriddedPart3d::generateGeometry( const std::array& input nVertCells += layersPerRegion[region].size(); } - const std::vector thicknessFactors = { -1.0, 0.0, 1.0 }; - const int nThicknessCells = 2; - cvf::Vec3d tVec = modelThickness * thicknessDirection; + const std::vector thicknessFactors = { -1.0, 0.0, 1.0 }; + const int nThicknessCells = 2; + cvf::Vec3d tVec = modelThickness * thicknessDirection; + std::vector> faultLines; + + for ( int i = 0; i <= nThicknessCells; i++ ) + { + faultLines.push_back( caf::Line( topBottomFaultPoints.first + thicknessFactors[i] * tVec, + topBottomFaultPoints.second + thicknessFactors[i] * tVec ) ); + } size_t reserveSize = ( nVertCells + 1 ) * ( nHorzCells + 1 ) * ( nThicknessCells + 1 ); m_nodes.reserve( reserveSize ); @@ -316,11 +327,26 @@ void RigGriddedPart3d::generateGeometry( const std::array& input for ( int h = 0; h <= (int)nHorzCells; h++ ) { - p = toPos - horizontalPartition[h] * stepHorz; + if ( h == (int)nHorzCells ) + { + p = toPos; + } + else + { + p = toPos - horizontalPartition[h] * stepHorz; + } for ( int t = 0; t <= nThicknessCells; t++, nodeIndex++ ) { - m_nodes.push_back( p + thicknessFactors[t] * tVec ); + auto nodePoint = p + thicknessFactors[t] * tVec; + + if ( ( h == (int)nHorzCells ) && + ( ( region == Regions::Reservoir ) || region == Regions::LowerOverburden || region == Regions::UpperUnderburden ) ) + { + nodePoint = faultLines[t].closestPointOnLine( nodePoint ); + } + + m_nodes.push_back( nodePoint ); if ( h == (int)nHorzCells ) { @@ -342,6 +368,11 @@ void RigGriddedPart3d::generateGeometry( const std::array& input else if ( h == (int)nHorzCells ) { m_boundaryNodes[Boundary::Fault].push_back( nodeIndex ); + + if ( region == Regions::Reservoir ) + { + m_boundaryNodes[Boundary::Reservoir].push_back( nodeIndex ); + } } } } @@ -767,3 +798,14 @@ void RigGriddedPart3d::updateElementSet( ElementSets elSet, } } } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigGriddedPart3d::postProcessBoundaryNodes( const caf::Line& faultLine ) +{ + for ( auto nodeIdx : m_boundaryNodes.at( RimFaultReactivation::Boundary::Reservoir ) ) + { + m_nodes[nodeIdx] = faultLine.closestPointOnLine( m_nodes[nodeIdx] ); + } +} diff --git a/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h b/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h index 49f44269c8..629e449d38 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h +++ b/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h @@ -24,6 +24,8 @@ #include "cvfObject.h" #include "cvfVector3.h" +#include "cafLine.h" + #include #include #include @@ -55,12 +57,14 @@ class RigGriddedPart3d : public cvf::Object double modelThickness, double topHeight, cvf::Vec3d thicknessDirection, - int nFaultZoneCells ); + int nFaultZoneCells, + std::pair topBottomFaultPoints ); void generateLocalNodes( const cvf::Mat4d transform ); void setUseLocalCoordinates( bool useLocalCoordinates ); void postProcessElementSets( const RigMainGrid* mainGrid, const RigActiveCellInfo* cellInfo ); + void postProcessBoundaryNodes( const caf::Line& faultLine ); bool useLocalCoordinates() const; double topHeight() const; diff --git a/Fwk/AppFwk/cafVizExtensions/cafLine.h b/Fwk/AppFwk/cafVizExtensions/cafLine.h index 2c549a2361..272c0548e5 100644 --- a/Fwk/AppFwk/cafVizExtensions/cafLine.h +++ b/Fwk/AppFwk/cafVizExtensions/cafLine.h @@ -56,6 +56,8 @@ class Line Line findLineBetweenNearestPoints( const Line& otherLine, bool* withinLineSegments = nullptr ); + cvf::Vector3 closestPointOnLine( const cvf::Vector3& point ) const; + private: cvf::Vector3 m_start; cvf::Vector3 m_end; diff --git a/Fwk/AppFwk/cafVizExtensions/cafLine.inl b/Fwk/AppFwk/cafVizExtensions/cafLine.inl index a0795abc07..02c110d431 100644 --- a/Fwk/AppFwk/cafVizExtensions/cafLine.inl +++ b/Fwk/AppFwk/cafVizExtensions/cafLine.inl @@ -105,3 +105,27 @@ caf::Line caf::Line::findLineBetweenNearestPoints( const Line& otherLine, } return Line( start() + s * d1, otherLine.start() + t * d2 ); } + +//-------------------------------------------------------------------------------------------------- +/// Returns the closest point on the line segment to the given point +/// From https://paulbourke.net/geometry/pointlineplane/ +//-------------------------------------------------------------------------------------------------- +template +cvf::Vector3 caf::Line::closestPointOnLine( const cvf::Vector3& point ) const +{ + S distance = vector().length(); + + const double u = ( ( ( point.x() - start().x() ) * ( end().x() - start().x() ) ) + + ( ( point.y() - start().y() ) * ( end().y() - start().y() ) ) + + ( ( point.z() - start().z() ) * ( end().z() - start().z() ) ) ) / + ( distance * distance ); + + if ( u < 0.0 ) return start(); + if ( u > 1.0 ) return end(); + + cvf::Vector3 projectedPoint( ( start().x() + u * ( end().x() - start().x() ) ), + ( start().y() + u * ( end().y() - start().y() ) ), + ( start().z() + u * ( end().z() - start().z() ) ) ); + + return projectedPoint; +}