Fault activation: node adjustments along fault (#11235)

Make sure same thickness vectors and fault lines are used for both parts.

ApplicationLibCode/FileInterface/RifFaultReactivationModelExporter.cpp

@@ -204,13 +204,13 @@ std::pair<bool, std::string> RifFaultReactivationModelExporter::printParts(
         for ( auto [boundary, boundaryName] : boundaries )
         {
             // Create boundary condition sets for each side of the parts (except top).
-            auto boundaryNodes    = grid->boundaryNodes();
+            const auto& boundaryNodes    = grid->boundaryNodes();
-            auto boundaryElements = grid->boundaryElements();
+            const auto& boundaryElements = grid->boundaryElements();
 
-            const std::vector<unsigned int>& nodes = boundaryNodes[boundary];
+            const std::vector<unsigned int>& nodes = boundaryNodes.at( boundary );
             RifInpExportTools::printNodeSet( stream, boundaryName, false, nodes );
 
-            const std::vector<unsigned int>& elements = boundaryElements[boundary];
+            const std::vector<unsigned int>& elements = boundaryElements.at( boundary );
             RifInpExportTools::printElementSet( stream, boundaryName, false, elements );
         }
 

ApplicationLibCode/ProjectDataModel/Faults/RimFaultReactivationEnums.h

@@ -39,7 +39,8 @@ enum class Boundary
 {
     FarSide,
     Bottom,
-    Fault
+    Fault,
+    Reservoir
 };
 
 enum class ElementSets

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 );
     }
 }

ApplicationLibCode/ReservoirDataModel/RigFaultReactivationModelGenerator.cpp

@@ -574,23 +574,38 @@ void RigFaultReactivationModelGenerator::generateGeometry( size_t            sta
 
     generatePointsFrontBack();
 
+    cvf::Vec3d                     tVec = m_modelThickness * m_modelNormal;
+    std::vector<cvf::Vec3d>        thicknessVectors;
+    std::vector<caf::Line<double>> faultLines;
+    const std::vector<double>      thicknessFactors = { -1.0, 0.0, 1.0 };
+
+    for ( int i = 0; i < 3; i++ )
+    {
+        faultLines.push_back( caf::Line<double>( 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 );

ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.cpp

@@ -25,8 +25,11 @@
 #include "RimFaultReactivationEnums.h"
 
 #include "cvfBoundingBox.h"
+#include "cvfPlane.h"
 #include "cvfTextureImage.h"
 
+#include "cafLine.h"
+
 #include <cmath>
 #include <map>
 
@@ -165,7 +168,7 @@ std::vector<double> RigGriddedPart3d::generateGrowingLayers( double zFrom, doubl
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-std::vector<double> RigGriddedPart3d::extractZValues( std::vector<cvf::Vec3d> points )
+std::vector<double> RigGriddedPart3d::extractZValues( const std::vector<cvf::Vec3d>& points )
 {
     std::vector<double> layers;
 
@@ -199,15 +202,15 @@ std::vector<double> RigGriddedPart3d::extractZValues( std::vector<cvf::Vec3d> po
 /// Assumes horizontal lines are parallel
 ///
 //--------------------------------------------------------------------------------------------------
-void RigGriddedPart3d::generateGeometry( const std::array<cvf::Vec3d, 12>& inputPoints,
+void RigGriddedPart3d::generateGeometry( const std::array<cvf::Vec3d, 12>&    inputPoints,
-                                         const std::vector<cvf::Vec3d>&    reservoirLayers,
+                                         const std::vector<cvf::Vec3d>&       reservoirLayers,
-                                         const double                      maxCellHeight,
+                                         double                               maxCellHeight,
-                                         double                            cellSizeFactor,
+                                         double                               cellSizeFactor,
-                                         const std::vector<double>&        horizontalPartition,
+                                         const std::vector<double>&           horizontalPartition,
-                                         double                            modelThickness,
+                                         const std::vector<caf::Line<double>> faultLines,
-                                         double                            topHeight,
+                                         const std::vector<cvf::Vec3d>&       thicknessVectors,
-                                         cvf::Vec3d                        thicknessDirection,
+                                         double                               topHeight,
-                                         int                               nFaultZoneCells )
+                                         int                                  nFaultZoneCells )
 {
     reset();
 
@@ -224,9 +227,10 @@ void RigGriddedPart3d::generateGeometry( const std::array<cvf::Vec3d, 12>& 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::Bottom]    = {};
-    m_boundaryNodes[Boundary::FarSide] = {};
+    m_boundaryNodes[Boundary::FarSide]   = {};
-    m_boundaryNodes[Boundary::Fault]   = {};
+    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<cvf::Vec3d, 12>& input
         nVertCells += layersPerRegion[region].size();
     }
 
-    const std::vector<double> thicknessFactors = { -1.0, 0.0, 1.0 };
+    const int nThicknessCells = 2;
-    const int                 nThicknessCells  = 2;
-    cvf::Vec3d                tVec             = modelThickness * thicknessDirection;
 
     size_t reserveSize = ( nVertCells + 1 ) * ( nHorzCells + 1 ) * ( nThicknessCells + 1 );
     m_nodes.reserve( reserveSize );
@@ -316,12 +318,31 @@ void RigGriddedPart3d::generateGeometry( const std::array<cvf::Vec3d, 12>& 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<cvf::Vec3d, 12>& input
                     else if ( h == (int)nHorzCells )
                     {
                         m_boundaryNodes[Boundary::Fault].push_back( nodeIndex );
+
+                        if ( region == Regions::Reservoir )
+                        {
+                            m_boundaryNodes[Boundary::Reservoir].push_back( nodeIndex );
+                        }
                     }
                 }
             }

ApplicationLibCode/ReservoirDataModel/RigGriddedPart3d.h

@@ -24,6 +24,8 @@
 #include "cvfObject.h"
 #include "cvfVector3.h"
 
+#include "cafLine.h"
+
 #include <array>
 #include <map>
 #include <set>
@@ -47,15 +49,15 @@ public:
 
     void reset();
 
-    void generateGeometry( const std::array<cvf::Vec3d, 12>& inputPoints,
+    void generateGeometry( const std::array<cvf::Vec3d, 12>&    inputPoints,
-                           const std::vector<cvf::Vec3d>&    reservoirLayers,
+                           const std::vector<cvf::Vec3d>&       reservoirLayers,
-                           double                            maxCellHeight,
+                           double                               maxCellHeight,
-                           double                            cellSizeFactor,
+                           double                               cellSizeFactor,
-                           const std::vector<double>&        horizontalPartition,
+                           const std::vector<double>&           horizontalPartition,
-                           double                            modelThickness,
+                           const std::vector<caf::Line<double>> faultLines,
-                           double                            topHeight,
+                           const std::vector<cvf::Vec3d>&       thicknessVectors,
-                           cvf::Vec3d                        thicknessDirection,
+                           double                               topHeight,
-                           int                               nFaultZoneCells );
+                           int                                  nFaultZoneCells );
 
     void generateLocalNodes( const cvf::Mat4d transform );
     void setUseLocalCoordinates( bool useLocalCoordinates );
@@ -87,7 +89,7 @@ protected:
     static cvf::Vec3d          stepVector( cvf::Vec3d start, cvf::Vec3d stop, int nSteps );
     static std::vector<double> generateConstantLayers( double zFrom, double zTo, double maxSize );
     static std::vector<double> generateGrowingLayers( double zFrom, double zTo, double maxSize, double growfactor );
-    static std::vector<double> extractZValues( std::vector<cvf::Vec3d> );
+    static std::vector<double> extractZValues( const std::vector<cvf::Vec3d>& points );
 
     void generateVerticalMeshlines( const std::vector<cvf::Vec3d>& cornerPoints, const std::vector<double>& horzPartition );
 

Fwk/AppFwk/cafVizExtensions/cafLine.h

@@ -56,6 +56,8 @@ public:
 
     Line findLineBetweenNearestPoints( const Line& otherLine, bool* withinLineSegments = nullptr );
 
+    cvf::Vector3<S> closestPointOnLine( const cvf::Vector3<S>& point ) const;
+
 private:
     cvf::Vector3<S> m_start;
     cvf::Vector3<S> m_end;

Fwk/AppFwk/cafVizExtensions/cafLine.inl

@@ -105,3 +105,27 @@ caf::Line<S> caf::Line<S>::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 <typename S>
+cvf::Vector3<S> caf::Line<S>::closestPointOnLine( const cvf::Vector3<S>& 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<S> projectedPoint( ( start().x() + u * ( end().x() - start().x() ) ),
+                                    ( start().y() + u * ( end().y() - start().y() ) ),
+                                    ( start().z() + u * ( end().z() - start().z() ) ) );
+
+    return projectedPoint;
+}