#9470 User defined perf length for azimuth fractures

Fractures: Add UI for user-defined perforation length for az fractures
Add user defined perforation length for azimuth fractures.
Janitor: remove duplication

ApplicationLibCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.cpp

@@ -46,6 +46,7 @@
 
 #include "RifEclipseSummaryAddress.h"
 #include "RifSummaryReaderInterface.h"
+
 #include "RigCaseCellResultsData.h"
 #include "RigEclipseCaseData.h"
 #include "RigEclipseToStimPlanCalculator.h"
@@ -62,6 +63,9 @@
 #include "RigWellPath.h"
 #include "RigWellPathStimplanIntersector.h"
 
+#include "cvfGeometryTools.h"
+#include "cvfPlane.h"
+
 #include <vector>
 
 //--------------------------------------------------------------------------------------------------
@@ -607,11 +611,13 @@ void RicExportFractureCompletionsImpl::calculateFractureToWellTransmissibilities
     gsl::not_null<const RigWellPath*>         wellPathGeometry,
     RigTransmissibilityCondenser&             transCondenser )
 {
-    ////
+    // If fracture has orientation Azimuth (without user-defined perforation length) or Transverse,
-    // If fracture has orientation Azimuth or Transverse, assume only radial inflow
+    // assume only radial inflow
+    bool useRadialInflow = ( fracTemplate->orientationType() == RimFractureTemplate::AZIMUTH &&
+                             !fracTemplate->useUserDefinedPerforationLength() ) ||
+                           fracTemplate->orientationType() == RimFractureTemplate::TRANSVERSE_WELL_PATH;
 
-    if ( fracTemplate->orientationType() == RimFractureTemplate::AZIMUTH ||
+    if ( useRadialInflow )
-         fracTemplate->orientationType() == RimFractureTemplate::TRANSVERSE_WELL_PATH )
     {
         std::pair<size_t, size_t> wellCellIJ = fractureGrid->fractureCellAtWellCenter();
         size_t wellCellIndex = fractureGrid->getGlobalIndexFromIJ( wellCellIJ.first, wellCellIJ.second );
@@ -630,12 +636,27 @@ void RicExportFractureCompletionsImpl::calculateFractureToWellTransmissibilities
                                                     { false, RigTransmissibilityCondenser::CellAddress::STIMPLAN, wellCellIndex },
                                                     radialTrans );
     }
-    else if ( fracTemplate->orientationType() == RimFractureTemplate::ALONG_WELL_PATH )
+    else
     {
-        ////
+        std::vector<cvf::Vec3d> wellPathPoints;
-        // If fracture has orientation along well, linear inflow along well and radial flow at endpoints
 
-        RigWellPathStimplanIntersector wellFractureIntersector( wellPathGeometry, fracture );
+        if ( fracTemplate->orientationType() == RimFractureTemplate::ALONG_WELL_PATH )
+        {
+            // If fracture has orientation along well, linear inflow along well and radial flow at endpoints
+            wellPathPoints =
+                wellPathGeometry->wellPathPointsIncludingInterpolatedIntersectionPoint( fracture->fractureMD() );
+        }
+        else
+        {
+            // If fracture has azimuth orientation with user-defined perforation length use linear inflow along
+            // well and radial flow at endpoints on a well path which is rotated into the fracture plane.
+            CAF_ASSERT( fracTemplate->orientationType() == RimFractureTemplate::AZIMUTH );
+            CAF_ASSERT( fracTemplate->useUserDefinedPerforationLength() );
+
+            wellPathPoints = computeWellPointsInFracturePlane( fracture, wellPathGeometry );
+        }
+
+        RigWellPathStimplanIntersector wellFractureIntersector( wellPathPoints, fracture );
         const std::map<size_t, RigWellPathStimplanIntersector::WellCellIntersection>& fractureWellCells =
             wellFractureIntersector.intersections();
 
@@ -671,6 +692,87 @@ void RicExportFractureCompletionsImpl::calculateFractureToWellTransmissibilities
     }
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+std::vector<cvf::Vec3d>
+    RicExportFractureCompletionsImpl::computeWellPointsInFracturePlane( gsl::not_null<const RimFracture*> fracture,
+                                                                        gsl::not_null<const RigWellPath*> wellPathGeometry )
+
+{
+    std::vector<cvf::Vec3d> wellPathPoints =
+        wellPathGeometry->wellPathPointsIncludingInterpolatedIntersectionPoint( fracture->fractureMD() );
+
+    RiaLogging::info( "Using user-defined perforation length on azimuth fracture." );
+    RiaLogging::info( QString( "Perforation length: %1" ).arg( fracture->perforationLength() ) );
+    double startMd = fracture->fractureMD() - fracture->perforationLength() / 2.0;
+    double endMd   = fracture->fractureMD() + fracture->perforationLength() / 2.0;
+
+    cvf::Vec3d anchorPosition = wellPathGeometry->interpolatedPointAlongWellPath( fracture->fractureMD() );
+
+    auto coordsAndMd = wellPathGeometry->clippedPointSubset( startMd, endMd );
+
+    auto       wellPathCoords = coordsAndMd.first;
+    cvf::Vec3d startPos       = wellPathCoords.front();
+    cvf::Vec3d endPos         = wellPathCoords.back();
+
+    cvf::Vec3d wellPathTangent = endPos - startPos;
+
+    cvf::Plane fracturePlane;
+    auto       fractureTransform = fracture->transformMatrix();
+    fracturePlane.setFromPointAndNormal( fractureTransform.translation(),
+                                         static_cast<cvf::Vec3d>( fractureTransform.col( 2 ) ) );
+
+    // Get the direction from the start position to the fracture plane
+    cvf::Vec3d startPosInFracturePlane            = fracturePlane.projectPoint( startPos );
+    cvf::Vec3d directionToStartPosInFracturePlane = ( startPosInFracturePlane - anchorPosition ).getNormalized();
+    cvf::Vec3d directionToStartPos                = ( startPos - anchorPosition ).getNormalized();
+
+    auto vecToString = []( const cvf::Vec3d& vec ) {
+        return QString( "[%1 %2 %3]" ).arg( vec.x() ).arg( vec.y() ).arg( vec.z() );
+    };
+
+    RiaLogging::info( QString( "Start perforation position: %1" ).arg( vecToString( startPos ) ) );
+    RiaLogging::info(
+        QString( "Start perforation position in fracture plane: %1" ).arg( vecToString( startPosInFracturePlane ) ) );
+    RiaLogging::info( QString( "Anchor pos: %1" ).arg( vecToString( anchorPosition ) ) );
+    RiaLogging::info(
+        QString( "Direction anchor to start position in plane: %1" ).arg( vecToString( directionToStartPosInFracturePlane ) ) );
+    RiaLogging::info(
+        QString( "Direction anchor to original start position: %1" ).arg( vecToString( directionToStartPos ) ) );
+
+    // Find the angle between the real direction (tangential to well path) and the fracture plane
+    double angle = cvf::GeometryTools::getAngle( directionToStartPosInFracturePlane, directionToStartPos );
+    RiaLogging::info( QString( "Angle: %1 degrees." ).arg( cvf::Math::toDegrees( angle ) ) );
+    auto rotMat =
+        cvf::GeometryTools::rotationMatrixBetweenVectors( directionToStartPos, directionToStartPosInFracturePlane );
+
+    auto rotatePoint = []( const cvf::Vec3d& point, const cvf::Vec3d& offset, auto rotMat ) {
+        cvf::Vec3d p = point - offset;
+        p.transformPoint( rotMat );
+        p += offset;
+        return p;
+    };
+
+    // Rotate the well path points into the plane
+    for ( auto& r : wellPathPoints )
+    {
+        r = rotatePoint( r, anchorPosition, rotMat );
+    }
+
+    cvf::Vec3d transformedStartPos = rotatePoint( startPos, anchorPosition, rotMat );
+    RiaLogging::info( QString( "Perforation start position: original: %1 --> adapted: %2" )
+                          .arg( vecToString( startPos ) )
+                          .arg( vecToString( transformedStartPos ) ) );
+
+    cvf::Vec3d transformedEndPos = rotatePoint( endPos, anchorPosition, rotMat );
+    RiaLogging::info( QString( "Perforation end position: original: %1 --> adapted: %2" )
+                          .arg( vecToString( endPos ) )
+                          .arg( vecToString( transformedEndPos ) ) );
+
+    return wellPathPoints;
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationLibCode/Commands/CompletionExportCommands/RicExportFractureCompletionsImpl.h

@@ -20,6 +20,8 @@
 
 #include "RigCompletionData.h"
 
+#include "cvfVector3.h"
+
 #include <gsl/gsl>
 
 #include <map>
@@ -151,5 +153,8 @@ private:
                                                const RigMainGrid*            mainGrid,
                                                const RigFractureGrid*        fractureGrid );
 
+    static std::vector<cvf::Vec3d> computeWellPointsInFracturePlane( gsl::not_null<const RimFracture*> fracture,
+                                                                     gsl::not_null<const RigWellPath*> wellPathGeometry );
+
     static bool loadResultsByName( RigCaseCellResultsData* cellResultsData, const std::vector<QString>& resultNames );
 };

ApplicationLibCode/ModelVisualization/RivWellConnectionFactorGeometryGenerator.cpp

@@ -21,6 +21,7 @@
 #include "cafEffectGenerator.h"
 #include "cvfArray.h"
 #include "cvfDrawableGeo.h"
+#include "cvfGeometryTools.h"
 #include "cvfPart.h"
 #include "cvfPrimitiveSetIndexedUInt.h"
 #include "cvfScalarMapper.h"
@@ -122,7 +123,8 @@ cvf::ref<cvf::DrawableGeo> RivWellConnectionFactorGeometryGenerator::createSurfa
 
     for ( const auto& item : m_completionVizData )
     {
-        auto rotMatrix = rotationMatrixBetweenVectors( cvf::Vec3d::Y_AXIS, item.m_direction );
+        auto rotMatrix =
+            cvf::Mat4f( cvf::GeometryTools::rotationMatrixBetweenVectors( cvf::Vec3d::Y_AXIS, item.m_direction ) );
 
         cvf::uint indexOffset = static_cast<cvf::uint>( vertices->size() );
 
@@ -180,28 +182,6 @@ size_t RivWellConnectionFactorGeometryGenerator::mapFromTriangleToConnectionInde
     return connectionIndex;
 }
 
-//--------------------------------------------------------------------------------------------------
-/// Taken from OverlayNavigationCube::computeNewUpVector
-/// Consider move to geometry util class
-//--------------------------------------------------------------------------------------------------
-cvf::Mat4f RivWellConnectionFactorGeometryGenerator::rotationMatrixBetweenVectors( const cvf::Vec3d& v1,
-                                                                                   const cvf::Vec3d& v2 )
-{
-    using namespace cvf;
-
-    Vec3d rotAxis = v1 ^ v2;
-    rotAxis.normalize();
-
-    // Guard acos against out-of-domain input
-    const double dotProduct = Math::clamp( v1 * v2, -1.0, 1.0 );
-    const double angle      = Math::acos( dotProduct );
-    Mat4d        rotMat     = Mat4d::fromRotation( rotAxis, angle );
-
-    Mat4f myMat( rotMat );
-
-    return myMat;
-}
-
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationLibCode/ModelVisualization/RivWellConnectionFactorGeometryGenerator.h

@@ -68,11 +68,10 @@ private:
     size_t                     mapFromTriangleToConnectionIndex( cvf::uint triangleIndex ) const;
     cvf::ref<cvf::DrawableGeo> createSurfaceGeometry();
 
-    static cvf::Mat4f rotationMatrixBetweenVectors( const cvf::Vec3d& v1, const cvf::Vec3d& v2 );
+    static void createStarGeometry( std::vector<cvf::Vec3f>* vertices,
-    static void       createStarGeometry( std::vector<cvf::Vec3f>* vertices,
+                                    std::vector<cvf::uint>*  indices,
-                                          std::vector<cvf::uint>*  indices,
+                                    float                    radius,
-                                          float                    radius,
+                                    float                    thickness );
-                                          float                    thickness );
 
     static void createSimplifiedStarGeometry( std::vector<cvf::Vec3f>* vertices,
                                               std::vector<cvf::uint>*  indices,

ApplicationLibCode/ProjectDataModel/Completions/RimEllipseFractureTemplate.cpp

@@ -489,6 +489,7 @@ void RimEllipseFractureTemplate::defineUiOrdering( QString uiConfigName, caf::Pd
         group->add( &m_permeability );
         group->add( &m_width );
         group->add( &m_skinFactor );
+        group->add( &m_userDefinedPerforationLength );
         group->add( &m_perforationLength );
         group->add( &m_perforationEfficiency );
         group->add( &m_wellDiameter );

ApplicationLibCode/ProjectDataModel/Completions/RimFractureTemplate.cpp

@@ -133,6 +133,10 @@ RimFractureTemplate::RimFractureTemplate()
                                  caf::AppEnum<FracOrientationEnum>( TRANSVERSE_WELL_PATH ),
                                  "Fracture Orientation" );
 
+    CAF_PDM_InitScriptableField( &m_userDefinedPerforationLength,
+                                 "UserDefinedPerforationLength",
+                                 false,
+                                 "User-defined Perforation Length" );
     CAF_PDM_InitScriptableField( &m_azimuthAngle, "AzimuthAngle", 0.0f, "Azimuth Angle" );
 
     CAF_PDM_InitField( &m_skinFactor, "SkinFactor", 0.0f, "Skin Factor" );
@@ -496,10 +500,12 @@ void RimFractureTemplate::prepareFieldsForUiDisplay()
          m_orientationType == RimFractureTemplate::TRANSVERSE_WELL_PATH )
     {
         m_azimuthAngle.uiCapability()->setUiHidden( true );
+        m_userDefinedPerforationLength.uiCapability()->setUiHidden( true );
     }
     else if ( m_orientationType == RimFractureTemplate::AZIMUTH )
     {
         m_azimuthAngle.uiCapability()->setUiHidden( false );
+        m_userDefinedPerforationLength.uiCapability()->setUiHidden( false );
     }
 
     if ( m_orientationType == RimFractureTemplate::ALONG_WELL_PATH )
@@ -510,7 +516,10 @@ void RimFractureTemplate::prepareFieldsForUiDisplay()
     else
     {
         m_perforationEfficiency.uiCapability()->setUiHidden( true );
-        m_perforationLength.uiCapability()->setUiHidden( true );
+
+        bool hidePerforationLength =
+            !( m_orientationType == RimFractureTemplate::AZIMUTH && m_userDefinedPerforationLength() );
+        m_perforationLength.uiCapability()->setUiHidden( hidePerforationLength );
     }
 
     if ( m_conductivityType == FINITE_CONDUCTIVITY )
@@ -939,6 +948,14 @@ double RimFractureTemplate::perforationLength() const
     return m_perforationLength;
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+bool RimFractureTemplate::useUserDefinedPerforationLength() const
+{
+    return m_userDefinedPerforationLength;
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationLibCode/ProjectDataModel/Completions/RimFractureTemplate.h

@@ -138,6 +138,7 @@ public:
     double                                      wellDiameter() const;
     FracConductivityEnum                        conductivityType() const;
     double                                      perforationLength() const;
+    bool                                        useUserDefinedPerforationLength() const;
 
     double                            wellPathDepthAtFracture() const;
     virtual std::pair<double, double> wellPathDepthAtFractureRange() const  = 0;
@@ -215,6 +216,7 @@ protected:
     caf::PdmField<caf::AppEnum<FracOrientationEnum>>           m_orientationType;
     caf::PdmField<float>                                       m_azimuthAngle;
     caf::PdmField<float>                                       m_skinFactor;
+    caf::PdmField<bool>                                        m_userDefinedPerforationLength;
     caf::PdmField<double>                                      m_perforationLength;
     caf::PdmField<double>                                      m_perforationEfficiency;
     caf::PdmField<double>                                      m_wellDiameter;

ApplicationLibCode/ProjectDataModel/Completions/RimMeshFractureTemplate.cpp

@@ -28,6 +28,7 @@
 #include "RigFractureCell.h"
 #include "RigFractureGrid.h"
 #include "RigTransmissibilityEquations.h"
+#include "RigWellPath.h"
 #include "RigWellPathStimplanIntersector.h"
 
 #include "RimEclipseView.h"
@@ -275,7 +276,10 @@ WellFractureIntersectionData RimMeshFractureTemplate::wellFractureIntersectionDa
                 RiaWeightedMeanCalculator<double>  conductivityCalc;
                 RiaWeightedGeometricMeanCalculator betaFactorCalc;
 
-                RigWellPathStimplanIntersector intersector( rimWellPath->wellPathGeometry(), fractureInstance );
+                std::vector<cvf::Vec3d> wellPathPoints =
+                    rimWellPath->wellPathGeometry()->wellPathPointsIncludingInterpolatedIntersectionPoint(
+                        fractureInstance->fractureMD() );
+                RigWellPathStimplanIntersector intersector( wellPathPoints, fractureInstance );
                 for ( const auto& v : intersector.intersections() )
                 {
                     size_t fractureGlobalCellIndex = v.first;
@@ -498,6 +502,7 @@ void RimMeshFractureTemplate::defineUiOrdering( QString uiConfigName, caf::PdmUi
         group->add( &m_conductivityResultNameOnFile );
         group->add( &m_conductivityType );
         group->add( &m_skinFactor );
+        group->add( &m_userDefinedPerforationLength );
         group->add( &m_perforationLength );
         group->add( &m_perforationEfficiency );
         group->add( &m_wellDiameter );

ApplicationLibCode/ReservoirDataModel/Completions/RigWellPathStimplanIntersector.cpp

@@ -36,11 +36,9 @@
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
-RigWellPathStimplanIntersector::RigWellPathStimplanIntersector( gsl::not_null<const RigWellPath*> wellPathGeom,
+RigWellPathStimplanIntersector::RigWellPathStimplanIntersector( const std::vector<cvf::Vec3d>&    wellPathPoints,
                                                                 gsl::not_null<const RimFracture*> rimFracture )
 {
-    std::vector<cvf::Vec3d> wellPathPoints =
-        wellPathGeom->wellPathPointsIncludingInterpolatedIntersectionPoint( rimFracture->fractureMD() );
     cvf::Mat4d fractureXf = rimFracture->transformMatrix();
     double     wellRadius = rimFracture->wellRadius();
 

ApplicationLibCode/ReservoirDataModel/Completions/RigWellPathStimplanIntersector.h

@@ -50,7 +50,7 @@ public:
         double computeLength() const { return cvf::Math::sqrt( hlength * hlength + vlength * vlength ); }
     };
 
-    RigWellPathStimplanIntersector( gsl::not_null<const RigWellPath*> wellpathGeom,
+    RigWellPathStimplanIntersector( const std::vector<cvf::Vec3d>&    wellPathPoints,
                                     gsl::not_null<const RimFracture*> rimFracture );
 
     const std::map<size_t, WellCellIntersection>& intersections() const;

ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.cpp

@@ -445,24 +445,6 @@ void RigThermalFractureResultUtil::appendDataToResultStatistics( std::shared_ptr
     }
 }
 
-//--------------------------------------------------------------------------------------------------
-/// Taken from OverlayNavigationCube::computeNewUpVector
-/// Consider move to geometry util class
-//--------------------------------------------------------------------------------------------------
-cvf::Mat4d RigThermalFractureResultUtil::rotationMatrixBetweenVectors( const cvf::Vec3d& v1, const cvf::Vec3d& v2 )
-{
-    using namespace cvf;
-
-    Vec3d rotAxis = v1 ^ v2;
-    rotAxis.normalize();
-
-    // Guard acos against out-of-domain input
-    const double dotProduct = Math::clamp( v1 * v2, -1.0, 1.0 );
-    const double angle      = Math::acos( dotProduct );
-    Mat4d        rotMat     = Mat4d::fromRotation( rotAxis, angle );
-    return rotMat;
-}
-
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------
@@ -493,7 +475,7 @@ std::vector<cvf::Vec3d>
     plane.setFromPoints( p0, p1, p2 );
 
     cvf::Vec3d planeNormal = plane.normal().getNormalized();
-    auto       rotMat      = rotationMatrixBetweenVectors( planeNormal, ( cvf::Vec3d::Z_AXIS ).getNormalized() );
+    auto rotMat = cvf::GeometryTools::rotationMatrixBetweenVectors( planeNormal, ( cvf::Vec3d::Z_AXIS ).getNormalized() );
 
     for ( auto& r : relativePos )
     {
@@ -533,7 +515,7 @@ std::vector<cvf::Vec3d>
     // Make sure normal is pointing down
     if ( directionNormal.y() > 0.0 ) directionNormal *= -1.0;
 
-    auto rotMat2 = rotationMatrixBetweenVectors( directionNormal, cvf::Vec3d::X_AXIS );
+    auto rotMat2 = cvf::GeometryTools::rotationMatrixBetweenVectors( directionNormal, cvf::Vec3d::X_AXIS );
     for ( auto& r : relativePos )
     {
         r.transformVector( rotMat2 );

ApplicationLibCode/ReservoirDataModel/RigThermalFractureResultUtil.h

@@ -96,8 +96,6 @@ private:
 
     static double linearSampling( double minValue, double maxValue, int numSamples, std::vector<double>& samples );
 
-    static cvf::Mat4d rotationMatrixBetweenVectors( const cvf::Vec3d& v1, const cvf::Vec3d& v2 );
-
     static std::vector<cvf::Vec3d>
         getRelativeCoordinates( std::shared_ptr<const RigThermalFractureDefinition> fractureDefinition,
                                 size_t                                              timeStepIndex );

ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.cpp

@@ -165,6 +165,21 @@ double GeometryTools::getAngle( const cvf::Vec3d& v1, const cvf::Vec3d& v2 )
     return angle;
 }
 
+//--------------------------------------------------------------------------------------------------
+///
+//--------------------------------------------------------------------------------------------------
+cvf::Mat4d GeometryTools::rotationMatrixBetweenVectors( const cvf::Vec3d& v1, const cvf::Vec3d& v2 )
+{
+    cvf::Vec3d rotAxis = v1 ^ v2;
+    rotAxis.normalize();
+
+    // Guard acos against out-of-domain input
+    const double dotProduct = cvf::Math::clamp( v1 * v2, -1.0, 1.0 );
+    const double angle      = cvf::Math::acos( dotProduct );
+    cvf::Mat4d   rotMat     = cvf::Mat4d::fromRotation( rotAxis, angle );
+    return rotMat;
+}
+
 //--------------------------------------------------------------------------------------------------
 ///
 //--------------------------------------------------------------------------------------------------

ApplicationLibCode/ReservoirDataModel/cvfGeometryTools.h

@@ -22,6 +22,7 @@
 #include "cvfArray.h"
 #include "cvfArrayWrapperConst.h"
 #include "cvfMatrix3.h"
+#include "cvfMatrix4.h"
 #include <list>
 #include <map>
 
@@ -140,6 +141,8 @@ public:
                                                const std::vector<bool>& faceOverlapPolygonWindingSameAsCubeFaceFlags,
                                                std::vector<IndexType>*  partialFacePolygon,
                                                bool*                    m_partiallyFreeCubeFaceHasHoles );
+
+    static cvf::Mat4d rotationMatrixBetweenVectors( const cvf::Vec3d& v1, const cvf::Vec3d& v2 );
 };
 
 template <typename IndexType>