From 86c33c1b460ce4a8c8d4b53a218934b150778a4b Mon Sep 17 00:00:00 2001 From: jonjenssen <69144954+jonjenssen@users.noreply.github.com> Date: Tue, 27 Feb 2024 17:05:48 +0100 Subject: [PATCH] Fault activation: node adjustments along fault (#11235) Make sure same thickness vectors and fault lines are used for both parts. --- .../RifFaultReactivationModelExporter.cpp | 8 +-- .../Faults/RimFaultReactivationEnums.h | 3 +- .../RigFaultReactivationModel.cpp | 3 +- .../RigFaultReactivationModelGenerator.cpp | 23 +++++-- .../ReservoirDataModel/RigGriddedPart3d.cpp | 62 +++++++++++++------ .../ReservoirDataModel/RigGriddedPart3d.h | 22 ++++--- Fwk/AppFwk/cafVizExtensions/cafLine.h | 2 + Fwk/AppFwk/cafVizExtensions/cafLine.inl | 24 +++++++ 8 files changed, 109 insertions(+), 38 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..24250b086e 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModel.cpp @@ -261,6 +261,7 @@ void RigFaultReactivationModel::postProcessElementSets( const RimEclipseCase* eC for ( auto part : allGridParts() ) { - m_3dparts[part]->postProcessElementSets( eCase->mainGrid(), cellInfo ); + auto gridPart = m_3dparts[part]; + gridPart->postProcessElementSets( eCase->mainGrid(), cellInfo ); } } diff --git a/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp b/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp index aac975721e..d21b9d7fd4 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp @@ -574,23 +574,38 @@ void RigFaultReactivationModelGenerator::generateGeometry( size_t sta generatePointsFrontBack(); + cvf::Vec3d tVec = m_modelThickness * m_modelNormal; + std::vector thicknessVectors; + std::vector> faultLines; + const std::vector thicknessFactors = { -1.0, 0.0, 1.0 }; + + for ( int i = 0; i < 3; i++ ) + { + faultLines.push_back( caf::Line( m_topFault + thicknessFactors[i] * tVec, m_bottomFault + thicknessFactors[i] * tVec ) ); + thicknessVectors.push_back( thicknessFactors[i] * tVec ); + } + frontPart->generateGeometry( m_frontPoints, frontReservoirLayers, m_maxCellHeight, m_cellSizeHeightFactor, m_horizontalPartition, - m_modelThickness, + faultLines, + thicknessVectors, m_topReservoirFront.z(), - m_modelNormal, m_faultZoneCells ); + + std::reverse( faultLines.begin(), faultLines.end() ); + std::reverse( thicknessVectors.begin(), thicknessVectors.end() ); + backPart->generateGeometry( m_backPoints, backReservoirLayers, m_maxCellHeight, m_cellSizeHeightFactor, m_horizontalPartition, - m_modelThickness, + faultLines, + thicknessVectors, m_topReservoirBack.z(), - -1.0 * m_modelNormal, m_faultZoneCells ); frontPart->generateLocalNodes( m_localCoordTransform ); diff --git a/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp b/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp index f98e3c19f8..43566a56eb 100644 --- a/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp +++ b/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp @@ -25,8 +25,11 @@ #include "RimFaultReactivationEnums.h" #include "cvfBoundingBox.h" +#include "cvfPlane.h" #include "cvfTextureImage.h" +#include "cafLine.h" + #include #include @@ -165,7 +168,7 @@ std::vector RigGriddedPart3d::generateGrowingLayers( double zFrom, doubl //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- -std::vector RigGriddedPart3d::extractZValues( std::vector points ) +std::vector RigGriddedPart3d::extractZValues( const std::vector& points ) { std::vector layers; @@ -199,15 +202,15 @@ std::vector RigGriddedPart3d::extractZValues( std::vector po /// Assumes horizontal lines are parallel /// //-------------------------------------------------------------------------------------------------- -void RigGriddedPart3d::generateGeometry( const std::array& inputPoints, - const std::vector& reservoirLayers, - const double maxCellHeight, - double cellSizeFactor, - const std::vector& horizontalPartition, - double modelThickness, - double topHeight, - cvf::Vec3d thicknessDirection, - int nFaultZoneCells ) +void RigGriddedPart3d::generateGeometry( const std::array& inputPoints, + const std::vector& reservoirLayers, + double maxCellHeight, + double cellSizeFactor, + const std::vector& horizontalPartition, + const std::vector> faultLines, + const std::vector& thicknessVectors, + double topHeight, + int nFaultZoneCells ) { 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,7 @@ 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 int nThicknessCells = 2; size_t reserveSize = ( nVertCells + 1 ) * ( nHorzCells + 1 ) * ( nThicknessCells + 1 ); m_nodes.reserve( reserveSize ); @@ -316,12 +318,31 @@ 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 + thicknessVectors[t]; + // adjust points along the fault line inside the reservoir to make sure they end up at the fault + if ( ( h == (int)nHorzCells ) && + ( ( region == Regions::Reservoir ) || region == Regions::LowerOverburden || region == Regions::UpperUnderburden ) ) + { + cvf::Plane zPlane; + zPlane.setFromPointAndNormal( nodePoint, cvf::Vec3d::Z_AXIS ); + zPlane.intersect( faultLines[t].start(), faultLines[t].end(), &nodePoint ); + } + + m_nodes.push_back( nodePoint ); + + // move nodes at fault used for data extraction a bit away from the fault if ( h == (int)nHorzCells ) { m_dataNodes.push_back( p + safetyOffset ); @@ -342,6 +363,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 ); + } } } } diff --git a/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h b/ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h index 49f44269c8..554cea6400 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 @@ -47,15 +49,15 @@ class RigGriddedPart3d : public cvf::Object void reset(); - void generateGeometry( const std::array& inputPoints, - const std::vector& reservoirLayers, - double maxCellHeight, - double cellSizeFactor, - const std::vector& horizontalPartition, - double modelThickness, - double topHeight, - cvf::Vec3d thicknessDirection, - int nFaultZoneCells ); + void generateGeometry( const std::array& inputPoints, + const std::vector& reservoirLayers, + double maxCellHeight, + double cellSizeFactor, + const std::vector& horizontalPartition, + const std::vector> faultLines, + const std::vector& thicknessVectors, + double topHeight, + int nFaultZoneCells ); void generateLocalNodes( const cvf::Mat4d transform ); void setUseLocalCoordinates( bool useLocalCoordinates ); @@ -87,7 +89,7 @@ class RigGriddedPart3d : public cvf::Object static cvf::Vec3d stepVector( cvf::Vec3d start, cvf::Vec3d stop, int nSteps ); static std::vector generateConstantLayers( double zFrom, double zTo, double maxSize ); static std::vector generateGrowingLayers( double zFrom, double zTo, double maxSize, double growfactor ); - static std::vector extractZValues( std::vector ); + static std::vector extractZValues( const std::vector& points ); void generateVerticalMeshlines( const std::vector& cornerPoints, const std::vector& horzPartition ); 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; +}