#3393 Continous sensible tangent calculation for targets wo tangent constraint.

Using abs(90 -inclination) and inverse distance as weights.
Moved well path calculations into a separate class

ApplicationCode/Application/Tools/CMakeLists_files.cmake

@@ -33,6 +33,7 @@ ${CMAKE_CURRENT_LIST_DIR}/RiaWellPlanCalculator.h
 ${CMAKE_CURRENT_LIST_DIR}/RiaSCurveCalculator.h
 ${CMAKE_CURRENT_LIST_DIR}/RiaArcCurveCalculator.h
 ${CMAKE_CURRENT_LIST_DIR}/RiaJCurveCalculator.h
+${CMAKE_CURRENT_LIST_DIR}/RiaLineArcWellPathCalculator.h
 ${CMAKE_CURRENT_LIST_DIR}/RiaOffshoreSphericalCoords.h
 ${CMAKE_CURRENT_LIST_DIR}/RiaWeightedMeanCalculator.h
 ${CMAKE_CURRENT_LIST_DIR}/RiaWeightedMeanCalculator.inl
@@ -75,6 +76,7 @@ ${CMAKE_CURRENT_LIST_DIR}/RiaWellPlanCalculator.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RiaSCurveCalculator.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RiaArcCurveCalculator.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RiaJCurveCalculator.cpp
+${CMAKE_CURRENT_LIST_DIR}/RiaLineArcWellPathCalculator.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RiaWeightedGeometricMeanCalculator.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RiaWeightedHarmonicMeanCalculator.cpp
 ${CMAKE_CURRENT_LIST_DIR}/RiaOptionItemFactory.cpp

ApplicationCode/Application/Tools/RiaLineArcWellPathCalculator.cpp

@@ -0,0 +1,281 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2018- 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 "RiaLineArcWellPathCalculator.h"
+#include "cvfBase.h"
+#include "cvfAssert.h"
+#include "RiaOffshoreSphericalCoords.h"
+#include "RiaJCurveCalculator.h"
+#include "RiaSCurveCalculator.h"
+
+#define M_PI       3.14159265358979323846   // pi
+
+cvf::Vec3d smootheningTargetTangent(const cvf::Vec3d& p1, const cvf::Vec3d& p2, const cvf::Vec3d& p3);
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+RiaLineArcWellPathCalculator::RiaLineArcWellPathCalculator(const cvf::Vec3d& referencePointXyz,
+                                                           const std::vector<WellTarget>& activeWellPathTargets)
+{
+    // Handle incomplete input
+
+    if (activeWellPathTargets.size() < 2)
+    {
+        m_startTangent = cvf::Vec3d::ZERO;
+ 
+        if (activeWellPathTargets.size() == 1)
+        {
+            m_lineArcEndpoints.push_back(  activeWellPathTargets[0].targetPointXYZ + referencePointXyz );
+            m_targetStatuses.resize(activeWellPathTargets.size(), 
+                                    { !activeWellPathTargets[0].isTangentConstrained, 0.0, 0.0,
+                                       true, std::numeric_limits<double>::infinity(),
+                                       true, std::numeric_limits<double>::infinity() });
+        }
+
+        return;
+    }
+
+    m_targetStatuses.resize(activeWellPathTargets.size(), { false, 0.0, 0.0, 
+                                                            false, std::numeric_limits<double>::infinity(), 
+                                                            false, std::numeric_limits<double>::infinity() });
+
+    std::vector<WellTarget> adjustedWellPathTargets = activeWellPathTargets;
+
+    // Calculate sensible tangents for targets without a fixed one
+
+    if ( activeWellPathTargets.size() > 2 )
+    {
+        for ( size_t tIdx = 0; tIdx < activeWellPathTargets.size() - 2; ++tIdx )
+        {
+            if ( !activeWellPathTargets[tIdx+1].isTangentConstrained )
+            {
+                cvf::Vec3d tangent = smootheningTargetTangent(activeWellPathTargets[tIdx  ].targetPointXYZ,
+                                                              activeWellPathTargets[tIdx+1].targetPointXYZ,
+                                                              activeWellPathTargets[tIdx+2].targetPointXYZ);
+                RiaOffshoreSphericalCoords tangentSphCS(tangent);
+                adjustedWellPathTargets[tIdx+1].azimuth = tangentSphCS.azi();
+                adjustedWellPathTargets[tIdx+1].inclination = tangentSphCS.inc();
+                adjustedWellPathTargets[tIdx+1].isTangentConstrained = true;
+
+                m_targetStatuses[tIdx+1].hasDerivedTangent = true;
+                m_targetStatuses[tIdx+1].resultAzimuth = tangentSphCS.azi();
+                m_targetStatuses[tIdx+1].resultInclination = tangentSphCS.inc();
+
+            }
+        }
+    }
+
+
+    m_lineArcEndpoints.push_back(  activeWellPathTargets[0].targetPointXYZ + referencePointXyz );
+
+    // Handle first segment if it is not an S-Curve
+
+    size_t startSSegmentIdx = 0;
+    size_t endSSegementIdx = activeWellPathTargets.size() - 1;
+
+    if (!adjustedWellPathTargets[0].isTangentConstrained)
+    {
+        startSSegmentIdx = 1; 
+
+        const WellTarget& target1 = adjustedWellPathTargets[0];
+        const WellTarget& target2 = adjustedWellPathTargets[1];
+        WellTargetStatus& target1Status = m_targetStatuses[0];
+        WellTargetStatus& target2Status = m_targetStatuses[1];
+
+        if (adjustedWellPathTargets[1].isTangentConstrained)
+        {
+            // Create an upside down J curve from target 2 back to 1
+
+            RiaJCurveCalculator jCurve(target2.targetPointXYZ,
+                                       target2.azimuth + M_PI,
+                                       M_PI - target2.inclination,
+                                       target2.radius1,
+                                       target1.targetPointXYZ);
+
+            if ( jCurve.curveStatus() == RiaJCurveCalculator::OK )
+            {
+                m_lineArcEndpoints.push_back(jCurve.firstArcEndpoint() + referencePointXyz);
+            }
+            else if ( jCurve.curveStatus() == RiaJCurveCalculator::FAILED_RADIUS_TOO_LARGE )
+            {
+                target2Status.hasOverriddenRadius1 = true;
+                target2Status.resultRadius1 = jCurve.radius();
+            }
+
+            m_lineArcEndpoints.push_back(target2.targetPointXYZ + referencePointXyz);
+
+            target1Status.hasDerivedTangent = true;
+            target1Status.resultAzimuth = jCurve.endAzimuth() + M_PI;
+            target1Status.resultInclination = M_PI - jCurve.endInclination();
+        }
+        else // The complete wellpath is a straight line from target 1 to 2
+        {
+            m_lineArcEndpoints.push_back(target2.targetPointXYZ + referencePointXyz );
+            cvf::Vec3d t12 = target2.targetPointXYZ - target1.targetPointXYZ;
+            RiaOffshoreSphericalCoords t12Sph(t12);
+
+            target1Status.hasDerivedTangent = true;
+            target1Status.resultAzimuth = t12Sph.azi();
+            target1Status.resultInclination = t12Sph.inc();
+
+            target2Status.hasDerivedTangent = true;
+            target2Status.resultAzimuth = t12Sph.azi();
+            target2Status.resultInclination = t12Sph.inc();
+        }
+
+        m_startTangent = RiaOffshoreSphericalCoords::unitVectorFromAziInc( target1Status.resultAzimuth, target1Status.resultInclination);
+    }
+    else
+    {
+        m_startTangent = RiaOffshoreSphericalCoords::unitVectorFromAziInc( activeWellPathTargets[0].azimuth, activeWellPathTargets[0].inclination);
+    }
+
+    if (!adjustedWellPathTargets.back().isTangentConstrained)
+    {
+        endSSegementIdx -= 1;
+    }
+
+    // Calculate S-curves
+
+    if ( activeWellPathTargets.size() > 1 )
+    {
+        for ( size_t tIdx = startSSegmentIdx; tIdx < endSSegementIdx; ++tIdx )
+        {
+            const WellTarget& target1 = adjustedWellPathTargets[tIdx];
+            const WellTarget& target2 = adjustedWellPathTargets[tIdx+1];
+            WellTargetStatus& target1Status = m_targetStatuses[tIdx];
+            WellTargetStatus& target2Status = m_targetStatuses[tIdx+1];
+
+            // Ignore targets in the same place
+            if ( (target1.targetPointXYZ - target2.targetPointXYZ).length() < 1e-6 ) continue;
+
+            if ( target1.isTangentConstrained
+                && target2.isTangentConstrained )
+            {
+                RiaSCurveCalculator sCurveCalc(target1.targetPointXYZ,
+                                               target1.azimuth,
+                                               target1.inclination,
+                                               target1.radius2,
+                                               target2.targetPointXYZ,
+                                               target2.azimuth,
+                                               target2.inclination,
+                                               target2.radius1);
+
+                if ( sCurveCalc.solveStatus() != RiaSCurveCalculator::CONVERGED )
+                {
+                    double p1p2Length = (target2.targetPointXYZ - target1.targetPointXYZ).length();
+                    sCurveCalc = RiaSCurveCalculator::fromTangentsAndLength(target1.targetPointXYZ,
+                                                                            target1.azimuth,
+                                                                            target1.inclination,
+                                                                            0.2*p1p2Length,
+                                                                            target2.targetPointXYZ,
+                                                                            target2.azimuth,
+                                                                            target2.inclination,
+                                                                            0.2*p1p2Length);
+
+                    //RiaLogging::warning("Using fall-back calculation of well path geometry between active target number: " + QString::number(tIdx+1) + " and " + QString::number(tIdx+2));
+
+                    target1Status.hasOverriddenRadius2 = true;
+                    target1Status.resultRadius2 = sCurveCalc.firstRadius();
+
+                    target2Status.hasOverriddenRadius1 = true;
+                    target2Status.resultRadius1 = sCurveCalc.secondRadius();
+                }
+
+                m_lineArcEndpoints.push_back(sCurveCalc.firstArcEndpoint() + referencePointXyz);
+                m_lineArcEndpoints.push_back(sCurveCalc.secondArcStartpoint() + referencePointXyz);
+                m_lineArcEndpoints.push_back(target2.targetPointXYZ + referencePointXyz);
+            }
+
+        }
+    }
+
+    // Handle last segment if (its not the same as the first) and it has not been handled as an S-Curve
+
+    if ( adjustedWellPathTargets.size() > 2 && endSSegementIdx < (adjustedWellPathTargets.size() - 1) )
+    {
+        size_t targetCount = adjustedWellPathTargets.size();
+        const WellTarget& target1 = adjustedWellPathTargets[targetCount-2];
+        const WellTarget& target2 = adjustedWellPathTargets[targetCount-1];
+        WellTargetStatus& target1Status = m_targetStatuses[targetCount-2];
+        WellTargetStatus& target2Status = m_targetStatuses[targetCount-1];
+
+        // Create an ordinary J curve
+
+        RiaJCurveCalculator jCurve(target1.targetPointXYZ,
+                                   target1.azimuth,
+                                   target1.inclination,
+                                   target1.radius2,
+                                   target2.targetPointXYZ);
+
+        if ( jCurve.curveStatus() == RiaJCurveCalculator::OK )
+        {
+            m_lineArcEndpoints.push_back(jCurve.firstArcEndpoint() + referencePointXyz);
+        }
+        else if ( jCurve.curveStatus() == RiaJCurveCalculator::FAILED_RADIUS_TOO_LARGE )
+        {
+            target1Status.hasOverriddenRadius2 = true;
+            target1Status.resultRadius2 = jCurve.radius();
+        }
+
+        m_lineArcEndpoints.push_back(target2.targetPointXYZ + referencePointXyz);
+
+        target2Status.hasDerivedTangent = true;
+        target2Status.resultAzimuth = jCurve.endAzimuth();
+        target2Status.resultInclination = jCurve.endInclination();
+    }
+
+}
+
+
+cvf::Vec3d smootheningTargetTangent(const cvf::Vec3d& p1, const cvf::Vec3d& p2, const cvf::Vec3d& p3)
+{
+    cvf::Vec3d t12 = p2 - p1;
+    cvf::Vec3d t23 = p3 - p2;
+
+    double length12 = t12.length();
+    double length23 = t23.length();
+
+    t12 /= length12; // Normalize
+    t23 /= length23; // Normalize
+
+    cvf::Vec3d t1t2Hor(t12);
+    t1t2Hor.z() = 0.0;
+    double t12HorLength = t1t2Hor.length();
+
+    cvf::Vec3d t23Hor(t23);
+    t23Hor.z() = 0.0;
+    double t23HorLength = t23Hor.length();
+
+    // Calculate weights as combo of inverse distance and horizontal component
+
+    double w12 = t12HorLength * 1.0/length12;
+    double w23 = t23HorLength * 1.0/length23;
+
+    // Weight the tangents 
+
+    t12 *= w12; // Weight
+    t23 *= w23; // Weight
+
+                // Sum and normalization of weights
+    cvf::Vec3d averageTangent = 1.0/(w12 + w23) * (t12 + t23);
+
+    return averageTangent;
+}
+

ApplicationCode/Application/Tools/RiaLineArcWellPathCalculator.h

@@ -0,0 +1,63 @@
+/////////////////////////////////////////////////////////////////////////////////
+//
+//  Copyright (C) 2018- 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.
+//
+/////////////////////////////////////////////////////////////////////////////////
+#pragma once
+
+#include "cvfBase.h"
+#include "cvfVector3.h"
+
+#include <vector>
+
+class RiaLineArcWellPathCalculator
+{
+public:
+    struct WellTarget
+    {
+        cvf::Vec3d targetPointXYZ;
+        bool isTangentConstrained;
+        double azimuth;
+        double inclination;
+
+        double radius1;
+        double radius2;
+    };
+
+    RiaLineArcWellPathCalculator(const cvf::Vec3d& referencePointXyz, 
+                                 const std::vector<RiaLineArcWellPathCalculator::WellTarget>& targets);
+
+    struct WellTargetStatus
+    {
+        bool hasDerivedTangent;
+        double resultAzimuth;
+        double resultInclination;
+
+        bool hasOverriddenRadius1;
+        double resultRadius1;
+        bool hasOverriddenRadius2;
+        double resultRadius2;
+    };
+
+    cvf::Vec3d                           startTangent() const { return m_startTangent; }
+    const std::vector<cvf::Vec3d>&       lineArcEndpoints() const { return m_lineArcEndpoints;} 
+    const std::vector<WellTargetStatus>& targetStatuses() const { return m_targetStatuses;}
+ 
+private:
+    cvf::Vec3d                    m_startTangent;
+    std::vector<cvf::Vec3d>       m_lineArcEndpoints;
+    std::vector<WellTargetStatus> m_targetStatuses;
+};
+

ApplicationCode/Commands/WellPathCommands/RicCreateWellTargetsPickEventHandler.cpp

@@ -83,8 +83,6 @@ bool RicCreateWellTargetsPickEventHandler::handlePickEvent(const Ric3DPickEvent&
 
         m_geometryToAddTargetsTo->insertTarget(nullptr, newTarget);
 
-        m_geometryToAddTargetsTo->addSmootheningTangentToNextToLastTargetIfSensible();
-
         m_geometryToAddTargetsTo->updateConnectedEditors();
         m_geometryToAddTargetsTo->updateWellPathVisualization();
 

ApplicationCode/ProjectDataModel/RimWellPathGeometryDef.cpp

@@ -128,9 +128,12 @@ cvf::ref<RigWellPath> RimWellPathGeometryDef::createWellPathGeometry()
 {
     cvf::ref<RigWellPath> wellPathGeometry = new RigWellPath; 
     
-    if (activeWellTargets().size() < 2) return wellPathGeometry;
+    RiaLineArcWellPathCalculator wellPathCalculator = lineArcWellPathCalculator();
+
+    if (wellPathCalculator.lineArcEndpoints().size() < 2) return wellPathGeometry;
+
+    RiaPolyArcLineSampler arcLineSampler(wellPathCalculator.startTangent(),  wellPathCalculator.lineArcEndpoints());
 
-    RiaPolyArcLineSampler arcLineSampler(startTangent(),  lineArcEndpoints());
 
     arcLineSampler.sampledPointsAndMDs(30,
                                          false,
@@ -144,7 +147,10 @@ cvf::ref<RigWellPath> RimWellPathGeometryDef::createWellPathGeometry()
 //--------------------------------------------------------------------------------------------------
 std::vector<RiaWellPlanCalculator::WellPlanSegment> RimWellPathGeometryDef::wellPlan() const
 {
-    RiaWellPlanCalculator wpCalc(startTangent(),  lineArcEndpoints());
+    RiaLineArcWellPathCalculator wellPathCalculator = lineArcWellPathCalculator();
+
+    RiaWellPlanCalculator wpCalc(wellPathCalculator.startTangent(),  wellPathCalculator.lineArcEndpoints());
+
     return wpCalc.wellPlan();
 }
 
@@ -375,6 +381,7 @@ std::vector<RimWellPathTarget*> RimWellPathGeometryDef::activeWellTargets() cons
     return active;
 }
 
+#if 0 // Kept for reference a bit longer
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
@@ -528,21 +535,47 @@ std::vector<cvf::Vec3d> RimWellPathGeometryDef::lineArcEndpoints() const
     return endPoints;
 }
 
+#endif
+
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-cvf::Vec3d RimWellPathGeometryDef::startTangent() const
+RiaLineArcWellPathCalculator RimWellPathGeometryDef::lineArcWellPathCalculator() const
 {
-    std::vector<RimWellPathTarget*> wellTargets = activeWellTargets();
+    std::vector<RimWellPathTarget*>  wellTargets = activeWellTargets();
 
-    if (!wellTargets.empty() && wellTargets[0]->targetType() == RimWellPathTarget::POINT_AND_TANGENT)
+    std::vector< RiaLineArcWellPathCalculator::WellTarget> targetDatas;
+
+    for (auto wellTarget : wellTargets)
     {
-        return wellTargets[0]->tangent();
+        targetDatas.push_back(wellTarget->wellTargetData());
     }
-    else
+
+    RiaLineArcWellPathCalculator wellPathCalculator(referencePointXyz(), targetDatas);
+    const std::vector<RiaLineArcWellPathCalculator::WellTargetStatus>& targetStatuses = wellPathCalculator.targetStatuses();
+
+    for ( size_t tIdx = 0 ; tIdx < wellTargets.size(); ++tIdx )
     {
-        return { 0, 0, -1 };
+        wellTargets[tIdx]->flagRadius1AsIncorrect(false, 0 );
+        wellTargets[tIdx]->flagRadius2AsIncorrect(false, 0 );
+
+        if ( targetStatuses[tIdx].hasDerivedTangent )
+        {
+            wellTargets[tIdx]->setDerivedTangent(targetStatuses[tIdx].resultAzimuth, targetStatuses[tIdx].resultInclination);
+        }
+
+        if ( targetStatuses[tIdx].hasOverriddenRadius1 )
+        {
+            wellTargets[tIdx]->flagRadius1AsIncorrect(true, targetStatuses[tIdx].resultRadius1);
+        }
+
+        if ( targetStatuses[tIdx].hasOverriddenRadius2 )
+        {
+            wellTargets[tIdx]->flagRadius2AsIncorrect(true, targetStatuses[tIdx].resultRadius2);
+        }
     }
+
+    return wellPathCalculator;
 }
 
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ProjectDataModel/RimWellPathGeometryDef.h

@@ -26,6 +26,7 @@
 #include "cafPdmPtrField.h"
 #include "cafPdmChildArrayField.h"
 #include "RiaWellPlanCalculator.h"
+#include "RiaLineArcWellPathCalculator.h"
 
 
 class RimWellPath;
@@ -88,8 +89,7 @@ private:
     void                            initAfterRead() override;
     QList<caf::PdmOptionItemInfo>   calculateValueOptions(const caf::PdmFieldHandle* fieldNeedingOptions, bool* useOptionsOnly) override;
 
-    std::vector<cvf::Vec3d>         lineArcEndpoints() const;
+    RiaLineArcWellPathCalculator    lineArcWellPathCalculator() const;
-    cvf::Vec3d                      startTangent() const;
 
 private:
     caf::PdmField<cvf::Vec3d>                        m_referencePointUtmXyd;

ApplicationCode/ProjectDataModel/RimWellPathTarget.cpp

@@ -89,6 +89,23 @@ void RimWellPathTarget::setDerivedTangent(double azimuth, double inclination)
     }
 }
 
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+RiaLineArcWellPathCalculator::WellTarget RimWellPathTarget::wellTargetData()
+{
+    RiaLineArcWellPathCalculator::WellTarget targetData;
+    
+    targetData.targetPointXYZ = targetPointXYZ();
+    targetData.isTangentConstrained = (targetType() == POINT_AND_TANGENT);
+    targetData.azimuth = azimuth();
+    targetData.inclination = inclination();
+    targetData.radius1 = radius1();
+    targetData.radius2 = radius2();
+
+    return targetData;
+}
+
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
@@ -159,6 +176,8 @@ double RimWellPathTarget::radius1() const
     // Degrees pr 10m
 
     // Degrees pr 30m
+    if (fabs(m_dogleg1) < 1e-6) return std::numeric_limits<double>::infinity();
+
     return 30.0/cvf::Math::toRadians(m_dogleg1);
 }
 
@@ -172,6 +191,9 @@ double RimWellPathTarget::radius2() const
     // Degrees pr 10m
 
     // Degrees pr 30m
+
+    if (fabs(m_dogleg2) < 1e-6) return std::numeric_limits<double>::infinity();
+
     return 30.0/cvf::Math::toRadians(m_dogleg2);
 }
 

ApplicationCode/ProjectDataModel/RimWellPathTarget.h

@@ -24,6 +24,7 @@
 #include "cafAppEnum.h"
 #include "cafPdmField.h"
 #include "cafPdmCoreVec3d.h"
+#include "RiaLineArcWellPathCalculator.h"
 
 class RimWellPathTarget : public caf::PdmObject
 {
@@ -38,6 +39,8 @@ public:
     void setAsPointXYZAndTangentTarget(const cvf::Vec3d& point, double azimuth, double inclination);
     void setDerivedTangent(double azimuth, double inclination);
 
+    RiaLineArcWellPathCalculator::WellTarget wellTargetData();
+
     enum TargetTypeEnum { POINT_AND_TANGENT, POINT };
     TargetTypeEnum targetType() const;
     cvf::Vec3d     targetPointXYZ() const;