ResInsight/ApplicationLibCode/ReservoirDataModel/Completions/RigWellPathStimplanIntersector.cpp

/////////////////////////////////////////////////////////////////////////////////
//
//  Copyright (C) 2017-     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 <http://www.gnu.org/licenses/gpl.html>
//  for more details.
//
/////////////////////////////////////////////////////////////////////////////////

#include "RigWellPathStimplanIntersector.h"

#include "RigCellGeometryTools.h"
#include "RigFractureCell.h"
#include "RigFractureGrid.h"
#include "RigWellPath.h"

#include "RimFracture.h"
#include "RimFractureTemplate.h"
#include "RimSimWellFracture.h"
#include "RimStimPlanFractureTemplate.h"

#include "cvfMath.h"
#include "cvfMatrix4.h"

#include <cmath>

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

    std::vector<std::vector<cvf::Vec3d>> fractureGridCellPolygons;
    {
        if ( rimFracture->fractureGrid() )
        {
            const std::vector<RigFractureCell>& stpCells = rimFracture->fractureGrid()->fractureCells();
            for ( const auto& stpCell : stpCells )
            {
                fractureGridCellPolygons.push_back( stpCell.getPolygon() );
            }
        }
    }
    double perforationLength = rimFracture->perforationLength();

    calculate( fractureXf,
               wellPathPoints,
               wellRadius,
               perforationLength,
               fractureGridCellPolygons,
               m_stimPlanCellIdxToIntersectionInfoMap );
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const std::map<size_t, RigWellPathStimplanIntersector::RigWellPathStimplanIntersector::WellCellIntersection>&
    RigWellPathStimplanIntersector::intersections() const
{
    return m_stimPlanCellIdxToIntersectionInfoMap;
}

//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigWellPathStimplanIntersector::calculate( const cvf::Mat4d&                           fractureXf,
                                                const std::vector<cvf::Vec3d>&              wellPathPointsDomainCoords,
                                                double                                      wellRadius,
                                                double                                      perforationLength,
                                                const std::vector<std::vector<cvf::Vec3d>>& fractureGridCellPolygons,
                                                std::map<size_t, WellCellIntersection>& m_stimPlanCellIdxToIntersectionInfoMap )
{
    cvf::Mat4d toFractureXf = fractureXf.getInverted();

    std::vector<cvf::Vec3d> perforationLengthBoundingBoxPolygon;
    {
        double cicleRadius           = perforationLength / 2;
        int    pointsInCirclePolygon = 20;

        for ( int i = 0; i < pointsInCirclePolygon; i++ )
        {
            double x = cicleRadius * cvf::Math::cos( i * ( 2 * cvf::PI_D / pointsInCirclePolygon ) );
            double y = cicleRadius * cvf::Math::sin( i * ( 2 * cvf::PI_D / pointsInCirclePolygon ) );
            perforationLengthBoundingBoxPolygon.push_back( cvf::Vec3d( x, y, 0 ) );
        }
    }

    // Convert well path to fracture template system

    std::vector<cvf::Vec3d> fractureRelativeWellPathPoints;
    for ( const auto& wellPPoint : wellPathPointsDomainCoords )
    {
        fractureRelativeWellPathPoints.push_back( wellPPoint.getTransformedPoint( toFractureXf ) );
    }

    // Clip well path to fracture domain

    std::vector<std::vector<cvf::Vec3d>> wellPathPartsWithinFracture =
        RigCellGeometryTools::clipPolylineByPolygon( fractureRelativeWellPathPoints,
                                                     perforationLengthBoundingBoxPolygon,
                                                     RigCellGeometryTools::INTERPOLATE_LINE_Z );

    // Remove the part of the well path that is more than well radius away from the fracture plane

    std::vector<std::vector<cvf::Vec3d>> intersectingWellPathParts;

    for ( const auto& part : wellPathPartsWithinFracture )
    {
        std::vector<cvf::Vec3d> currentIntersectingWpPart;
        for ( size_t vxIdx = 0; vxIdx < part.size() - 1; ++vxIdx )
        {
            double thisAbsZ = fabs( part[vxIdx].z() );
            double nextAbsZ = fabs( part[vxIdx + 1].z() );
            double thisZ    = part[vxIdx].z();
            double nextZ    = part[vxIdx + 1].z();

            if ( thisAbsZ >= wellRadius && nextAbsZ >= wellRadius )
            {
                if ( ( thisZ >= 0 && nextZ >= 0 ) || ( thisZ <= 0 && nextZ <= 0 ) )
                {
                    continue; // Outside
                }
                else // In and out
                {
                    {
                        double wellRadiusDistFromPlane = thisZ > 0 ? wellRadius : -wellRadius;

                        double fraction = ( wellRadiusDistFromPlane - thisZ ) / ( nextZ - thisZ );

                        cvf::Vec3d intersectPoint = part[vxIdx] + fraction * ( part[vxIdx + 1] - part[vxIdx] );
                        currentIntersectingWpPart.push_back( intersectPoint );
                    }
                    {
                        double wellRadiusDistFromPlane = nextZ > 0 ? wellRadius : -wellRadius;

                        double fraction = ( wellRadiusDistFromPlane - thisZ ) / ( nextZ - thisZ );

                        cvf::Vec3d intersectPoint = part[vxIdx] + fraction * ( part[vxIdx + 1] - part[vxIdx] );
                        currentIntersectingWpPart.push_back( intersectPoint );

                        intersectingWellPathParts.push_back( currentIntersectingWpPart );
                        currentIntersectingWpPart.clear();
                    }
                    continue;
                }
            }
            if ( thisAbsZ < wellRadius && nextAbsZ < wellRadius ) // Inside
            {
                currentIntersectingWpPart.push_back( part[vxIdx] );
                continue;
            }

            if ( thisAbsZ < wellRadius && nextAbsZ >= wellRadius ) // Going out
            {
                currentIntersectingWpPart.push_back( part[vxIdx] );

                double wellRadiusDistFromPlane = nextZ > 0 ? wellRadius : -wellRadius;

                double fraction = ( wellRadiusDistFromPlane - thisZ ) / ( nextZ - thisZ );

                cvf::Vec3d intersectPoint = part[vxIdx] + fraction * ( part[vxIdx + 1] - part[vxIdx] );
                currentIntersectingWpPart.push_back( intersectPoint );

                intersectingWellPathParts.push_back( currentIntersectingWpPart );
                currentIntersectingWpPart.clear();
                continue;
            }

            if ( thisAbsZ >= wellRadius && nextAbsZ < wellRadius ) // Going in
            {
                double wellRadiusDistFromPlane = thisZ > 0 ? wellRadius : -wellRadius;

                double fraction = ( wellRadiusDistFromPlane - thisZ ) / ( nextZ - thisZ );

                cvf::Vec3d intersectPoint = part[vxIdx] + fraction * ( part[vxIdx + 1] - part[vxIdx] );
                currentIntersectingWpPart.push_back( intersectPoint );
                continue;
            }
        }

        // Add last point if it is within the radius

        if ( part.size() > 1 && fabs( part.back().z() ) < wellRadius )
        {
            currentIntersectingWpPart.push_back( part.back() );
        }

        if ( !currentIntersectingWpPart.empty() )
        {
            intersectingWellPathParts.push_back( currentIntersectingWpPart );
        }
    }

    // Find the StimPlan cells touched by the intersecting well path parts

    for ( size_t cIdx = 0; cIdx < fractureGridCellPolygons.size(); ++cIdx )
    {
        const std::vector<cvf::Vec3d>& cellPolygon = fractureGridCellPolygons[cIdx];
        for ( const auto& wellpathPart : intersectingWellPathParts )
        {
            std::vector<std::vector<cvf::Vec3d>> wellPathPartsInPolygon =
                RigCellGeometryTools::clipPolylineByPolygon( wellpathPart, cellPolygon, RigCellGeometryTools::USE_HUGEVAL );
            for ( const auto& wellPathPartInCell : wellPathPartsInPolygon )
            {
                if ( !wellPathPartInCell.empty() )
                {
                    int endpointCount = 0;
                    if ( wellPathPartInCell.front().z() != HUGE_VAL ) ++endpointCount;
                    if ( wellPathPartInCell.back().z() != HUGE_VAL ) ++endpointCount;

                    cvf::Vec3d intersectionLength = ( wellPathPartInCell.back() - wellPathPartInCell.front() );
                    double     xLengthInCell      = fabs( intersectionLength.x() );
                    double     yLengthInCell      = fabs( intersectionLength.y() );

                    m_stimPlanCellIdxToIntersectionInfoMap[cIdx].endpointCount += endpointCount;
                    m_stimPlanCellIdxToIntersectionInfoMap[cIdx].hlength += xLengthInCell;
                    m_stimPlanCellIdxToIntersectionInfoMap[cIdx].vlength += yLengthInCell;
                }
            }
        }
    }
}