mirror of
https://github.com/OPM/ResInsight.git
synced 2025-02-25 18:55:39 -06:00
Rename ApplicationCode to ApplicationLibCode
This commit is contained in:
Gaute Lindkvist
@@ -0,0 +1,34 @@
|
||||
|
||||
set (SOURCE_GROUP_HEADER_FILES
|
||||
|
||||
${CMAKE_CURRENT_LIST_DIR}/RiaPolyArcLineSampler.h
|
||||
${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
|
||||
)
|
||||
|
||||
set (SOURCE_GROUP_SOURCE_FILES
|
||||
${CMAKE_CURRENT_LIST_DIR}/RiaPolyArcLineSampler.cpp
|
||||
${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
|
||||
)
|
||||
|
||||
list(APPEND CODE_HEADER_FILES
|
||||
${SOURCE_GROUP_HEADER_FILES}
|
||||
)
|
||||
|
||||
list(APPEND CODE_SOURCE_FILES
|
||||
${SOURCE_GROUP_SOURCE_FILES}
|
||||
)
|
||||
|
||||
set (QT_MOC_HEADERS
|
||||
${QT_MOC_HEADERS}
|
||||
)
|
||||
|
||||
|
||||
source_group( "Application\\Tools\\WellPathTools" FILES ${SOURCE_GROUP_HEADER_FILES} ${SOURCE_GROUP_SOURCE_FILES} ${CMAKE_CURRENT_LIST_DIR}/CMakeLists_files.cmake )
|
||||
@@ -0,0 +1,102 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2018- Equinor 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 "RiaArcCurveCalculator.h"
|
||||
|
||||
#include "RiaOffshoreSphericalCoords.h"
|
||||
#include "cvfGeometryTools.h"
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// + p1
|
||||
/// t1 **
|
||||
/// * + C
|
||||
/// *
|
||||
/// + p2
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
RiaArcCurveCalculator::RiaArcCurveCalculator( cvf::Vec3d p1, cvf::Vec3d t1, cvf::Vec3d p2 )
|
||||
: m_radius( std::numeric_limits<double>::infinity() )
|
||||
, m_arcCS( cvf::Mat4d::ZERO )
|
||||
, m_endAzi( 0 )
|
||||
, m_endInc( 0 )
|
||||
, m_curveStatus( OK )
|
||||
{
|
||||
bool isOk = t1.normalize();
|
||||
if ( !isOk )
|
||||
{
|
||||
// No tangent. Bail out
|
||||
m_curveStatus = FAILED_INPUT_OVERLAP;
|
||||
return;
|
||||
}
|
||||
|
||||
cvf::Vec3d p1p2 = p2 - p1;
|
||||
cvf::Vec3d t12 = p1p2.getNormalized( &isOk );
|
||||
if ( !isOk )
|
||||
{
|
||||
// p1 and p2 in the same place.
|
||||
m_curveStatus = FAILED_INPUT_OVERLAP;
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
cvf::Vec3d N = ( t1 ^ t12 ).getNormalized( &isOk );
|
||||
if ( !isOk )
|
||||
{
|
||||
// P2 is on the p1 + k*t1 line. We have a straight line
|
||||
m_curveStatus = OK_STRAIGHT_LINE;
|
||||
|
||||
RiaOffshoreSphericalCoords endTangent( t1 );
|
||||
m_endTangent = t1;
|
||||
m_endAzi = endTangent.azi();
|
||||
m_endInc = endTangent.inc();
|
||||
m_radius = std::numeric_limits<double>::infinity();
|
||||
m_arcAngle = 0;
|
||||
m_arcLength = p1p2.length();
|
||||
return;
|
||||
}
|
||||
|
||||
cvf::Vec3d tr1 = ( N ^ t1 ).getNormalized();
|
||||
|
||||
m_radius = 0.5 * p1p2.length() / ( tr1.dot( t12 ) );
|
||||
|
||||
cvf::Vec3d C = p1 + m_radius * tr1;
|
||||
|
||||
cvf::Vec3d nTr1 = -tr1;
|
||||
m_arcCS = cvf::Mat4d::fromCoordSystemAxes( &nTr1, &t1, &N );
|
||||
|
||||
m_arcCS.setTranslation( C );
|
||||
|
||||
m_arcAngle = cvf::GeometryTools::getAngle( N, p1 - C, p2 - C );
|
||||
|
||||
m_arcLength = m_radius * m_arcAngle;
|
||||
|
||||
m_endTangent = N ^ ( p2 - C ).getNormalized();
|
||||
|
||||
RiaOffshoreSphericalCoords endTangent( m_endTangent );
|
||||
m_endAzi = endTangent.azi();
|
||||
m_endInc = endTangent.inc();
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
RiaArcCurveCalculator::RiaArcCurveCalculator( cvf::Vec3d p1, double azi1, double inc1, cvf::Vec3d p2 )
|
||||
{
|
||||
cvf::Vec3d t1( RiaOffshoreSphericalCoords::unitVectorFromAziInc( azi1, inc1 ) );
|
||||
|
||||
( *this ) = RiaArcCurveCalculator( p1, t1, p2 );
|
||||
}
|
||||
@@ -0,0 +1,68 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2018- Equinor 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 "cvfMatrix4.h"
|
||||
#include "cvfVector3.h"
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// + p1
|
||||
/// t1 **
|
||||
/// * + C
|
||||
/// *
|
||||
/// + p2
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
class RiaArcCurveCalculator
|
||||
{
|
||||
public:
|
||||
RiaArcCurveCalculator( cvf::Vec3d p1, cvf::Vec3d t1, cvf::Vec3d p2 );
|
||||
RiaArcCurveCalculator( cvf::Vec3d p1, double azi1, double inc1, cvf::Vec3d p2 );
|
||||
|
||||
enum CurveStatus
|
||||
{
|
||||
OK,
|
||||
OK_STRAIGHT_LINE,
|
||||
FAILED_INPUT_OVERLAP
|
||||
};
|
||||
|
||||
CurveStatus curveStatus() const { return m_curveStatus; }
|
||||
|
||||
cvf::Mat4d arcCS() const { return m_arcCS; }
|
||||
double radius() const { return m_radius; }
|
||||
double arcAngle() const { return m_arcAngle; }
|
||||
double arcLength() const { return m_arcLength; }
|
||||
cvf::Vec3d center() const { return m_arcCS.translation(); }
|
||||
cvf::Vec3d normal() const { return cvf::Vec3d( m_arcCS.col( 2 ) ); }
|
||||
|
||||
double endAzimuth() const { return m_endAzi; }
|
||||
double endInclination() const { return m_endInc; }
|
||||
cvf::Vec3d endTangent() const { return m_endTangent; }
|
||||
|
||||
private:
|
||||
CurveStatus m_curveStatus;
|
||||
|
||||
double m_radius;
|
||||
double m_arcLength;
|
||||
double m_arcAngle;
|
||||
cvf::Mat4d m_arcCS;
|
||||
|
||||
double m_endAzi;
|
||||
double m_endInc;
|
||||
cvf::Vec3d m_endTangent;
|
||||
};
|
||||
@@ -0,0 +1,95 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2018- Equinor 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 "RiaJCurveCalculator.h"
|
||||
#include "RiaArcCurveCalculator.h"
|
||||
#include "RiaOffshoreSphericalCoords.h"
|
||||
#include "cvfMatrix3.h"
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
RiaJCurveCalculator::RiaJCurveCalculator( cvf::Vec3d p1, double azi1, double inc1, double r1, cvf::Vec3d p2 )
|
||||
: m_c1( cvf::Vec3d::UNDEFINED )
|
||||
, m_n1( cvf::Vec3d::UNDEFINED )
|
||||
, m_radius( std::numeric_limits<double>::infinity() )
|
||||
, m_curveStatus( OK )
|
||||
{
|
||||
cvf::Vec3d t1( RiaOffshoreSphericalCoords::unitVectorFromAziInc( azi1, inc1 ) );
|
||||
|
||||
cvf::Vec3d p1p2 = p2 - p1;
|
||||
|
||||
cvf::Vec3d tr1 = p1p2 - ( p1p2.dot( t1 ) ) * t1;
|
||||
double tr1Length = tr1.length();
|
||||
if ( tr1Length < 1e-9 )
|
||||
{
|
||||
// p2 is on the p1 + t12 line. Degenerates to a line.
|
||||
m_curveStatus = OK_STRAIGHT_LINE;
|
||||
m_firstArcEndpoint = p2;
|
||||
m_endAzi = azi1;
|
||||
m_endInc = inc1;
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
tr1 /= tr1Length;
|
||||
|
||||
cvf::Vec3d c1 = p1 + r1 * tr1;
|
||||
cvf::Vec3d p2c1 = c1 - p2;
|
||||
|
||||
double p2c1Length = p2c1.length();
|
||||
if ( p2c1Length < r1 || r1 == std::numeric_limits<double>::infinity() )
|
||||
{
|
||||
// Radius is too big. We can not get to point 2 using the requested radius.
|
||||
m_curveStatus = FAILED_RADIUS_TOO_LARGE;
|
||||
|
||||
RiaArcCurveCalculator arc( p1, t1, p2 );
|
||||
if ( arc.curveStatus() == RiaArcCurveCalculator::OK || arc.curveStatus() == RiaArcCurveCalculator::OK_STRAIGHT_LINE )
|
||||
{
|
||||
m_c1 = arc.center();
|
||||
m_n1 = arc.normal();
|
||||
m_firstArcEndpoint = p2;
|
||||
m_endAzi = arc.endAzimuth();
|
||||
m_endInc = arc.endInclination();
|
||||
m_radius = arc.radius();
|
||||
}
|
||||
else
|
||||
{
|
||||
m_firstArcEndpoint = p2;
|
||||
m_endAzi = azi1;
|
||||
m_endInc = inc1;
|
||||
}
|
||||
return;
|
||||
}
|
||||
|
||||
double d = sqrt( p2c1Length * p2c1Length - r1 * r1 );
|
||||
|
||||
double betha = asin( r1 / p2c1Length );
|
||||
cvf::Vec3d tp2c1 = p2c1 / p2c1Length;
|
||||
cvf::Vec3d nc1 = t1 ^ tr1;
|
||||
|
||||
cvf::Vec3d tp11p2 = -tp2c1.getTransformedVector( cvf::Mat3d::fromRotation( nc1, betha ) );
|
||||
|
||||
m_firstArcEndpoint = p2 - d * tp11p2;
|
||||
m_c1 = c1;
|
||||
m_n1 = nc1;
|
||||
|
||||
RiaOffshoreSphericalCoords endTangent( tp11p2 );
|
||||
m_endAzi = endTangent.azi();
|
||||
m_endInc = endTangent.inc();
|
||||
}
|
||||
@@ -0,0 +1,67 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2018- Equinor 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 "cvfVector3.h"
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// + p1
|
||||
/// t1 **
|
||||
/// * r1 + C
|
||||
/// *
|
||||
/// + firstArcEndpoint
|
||||
/// *
|
||||
/// *
|
||||
/// + p2
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
class RiaJCurveCalculator
|
||||
{
|
||||
public:
|
||||
RiaJCurveCalculator( cvf::Vec3d p1, double azi1, double inc1, double r1, cvf::Vec3d p2 );
|
||||
enum CurveStatus
|
||||
{
|
||||
OK,
|
||||
OK_STRAIGHT_LINE,
|
||||
FAILED_INPUT_OVERLAP,
|
||||
FAILED_RADIUS_TOO_LARGE
|
||||
};
|
||||
|
||||
CurveStatus curveStatus() const { return m_curveStatus; }
|
||||
|
||||
cvf::Vec3d firstArcEndpoint() const { return m_firstArcEndpoint; }
|
||||
|
||||
double radius() const { return m_radius; }
|
||||
cvf::Vec3d firstCenter() const { return m_c1; }
|
||||
cvf::Vec3d firstNormal() const { return m_n1; }
|
||||
|
||||
double endAzimuth() const { return m_endAzi; }
|
||||
double endInclination() const { return m_endInc; }
|
||||
|
||||
private:
|
||||
CurveStatus m_curveStatus;
|
||||
|
||||
cvf::Vec3d m_firstArcEndpoint;
|
||||
|
||||
double m_radius;
|
||||
cvf::Vec3d m_c1;
|
||||
cvf::Vec3d m_n1;
|
||||
|
||||
double m_endAzi;
|
||||
double m_endInc;
|
||||
};
|
||||
@@ -0,0 +1,297 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2018- Equinor 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 "RiaJCurveCalculator.h"
|
||||
#include "RiaOffshoreSphericalCoords.h"
|
||||
#include "RiaSCurveCalculator.h"
|
||||
#include "cvfAssert.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,
|
||||
false,
|
||||
true,
|
||||
std::numeric_limits<double>::infinity(),
|
||||
false,
|
||||
true,
|
||||
std::numeric_limits<double>::infinity() } );
|
||||
}
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
m_targetStatuses.resize( activeWellPathTargets.size(),
|
||||
{ false,
|
||||
0.0,
|
||||
0.0,
|
||||
false,
|
||||
false,
|
||||
std::numeric_limits<double>::infinity(),
|
||||
false,
|
||||
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();
|
||||
|
||||
target2Status.isRadius1Editable = true;
|
||||
}
|
||||
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 );
|
||||
|
||||
target1Status.isRadius2Editable = true;
|
||||
target2Status.isRadius1Editable = true;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// 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 );
|
||||
|
||||
target1Status.isRadius2Editable = true;
|
||||
|
||||
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;
|
||||
}
|
||||
@@ -0,0 +1,65 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2018- Equinor 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 "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 isRadius1Editable;
|
||||
bool hasOverriddenRadius1;
|
||||
double resultRadius1;
|
||||
|
||||
bool isRadius2Editable;
|
||||
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;
|
||||
};
|
||||
@@ -0,0 +1,188 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2018- Equinor 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 "RiaPolyArcLineSampler.h"
|
||||
#include "RiaArcCurveCalculator.h"
|
||||
|
||||
#include "cvfGeometryTools.h"
|
||||
#include "cvfMatrix4.h"
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
RiaPolyArcLineSampler::RiaPolyArcLineSampler( const cvf::Vec3d& startTangent,
|
||||
const std::vector<cvf::Vec3d>& lineArcEndPoints )
|
||||
: m_startTangent( startTangent )
|
||||
, m_lineArcEndPoints( lineArcEndPoints )
|
||||
, m_maxSamplingsInterval( 0.15 )
|
||||
, m_isResamplingLines( true )
|
||||
, m_totalMD( 0.0 )
|
||||
{
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::pair<std::vector<cvf::Vec3d>, std::vector<double>>
|
||||
RiaPolyArcLineSampler::sampledPointsAndMDs( double sampleInterval, bool isResamplingLines )
|
||||
{
|
||||
CVF_ASSERT( sampleInterval > 0.0 );
|
||||
|
||||
m_maxSamplingsInterval = sampleInterval;
|
||||
m_isResamplingLines = isResamplingLines;
|
||||
|
||||
m_points.clear();
|
||||
m_meshDs.clear();
|
||||
|
||||
double startMD = 0.0;
|
||||
|
||||
std::vector<cvf::Vec3d> pointsNoDuplicates = RiaPolyArcLineSampler::pointsWithoutDuplicates( m_lineArcEndPoints );
|
||||
|
||||
if ( pointsNoDuplicates.size() < 2 ) return std::make_pair( m_points, m_meshDs );
|
||||
|
||||
m_totalMD = startMD;
|
||||
|
||||
cvf::Vec3d p1 = pointsNoDuplicates[0];
|
||||
|
||||
m_points.push_back( p1 );
|
||||
m_meshDs.push_back( m_totalMD );
|
||||
|
||||
cvf::Vec3d t2 = m_startTangent;
|
||||
|
||||
for ( size_t pIdx = 0; pIdx < pointsNoDuplicates.size() - 1; ++pIdx )
|
||||
{
|
||||
sampleSegment( t2, pointsNoDuplicates[pIdx], pointsNoDuplicates[pIdx + 1], &t2 );
|
||||
}
|
||||
|
||||
return std::make_pair( m_points, m_meshDs );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RiaPolyArcLineSampler::sampleSegment( cvf::Vec3d t1, cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent )
|
||||
{
|
||||
cvf::Vec3d p1p2 = p2 - p1;
|
||||
|
||||
CVF_ASSERT( p1p2.lengthSquared() > 1e-20 );
|
||||
|
||||
if ( cvf::GeometryTools::getAngle( t1, p1p2 ) < 1e-5 )
|
||||
{
|
||||
sampleLine( p1, p2, endTangent );
|
||||
}
|
||||
else // resample arc
|
||||
{
|
||||
sampleArc( t1, p1, p2, endTangent );
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::vector<cvf::Vec3d> RiaPolyArcLineSampler::pointsWithoutDuplicates( const std::vector<cvf::Vec3d>& points )
|
||||
{
|
||||
std::vector<cvf::Vec3d> outputPoints;
|
||||
|
||||
cvf::Vec3d previousPoint = cvf::Vec3d::UNDEFINED;
|
||||
const double threshold = 1e-6;
|
||||
for ( const auto& p : points )
|
||||
{
|
||||
if ( previousPoint.isUndefined() || ( ( previousPoint - p ).lengthSquared() ) > threshold )
|
||||
{
|
||||
outputPoints.push_back( p );
|
||||
previousPoint = p;
|
||||
}
|
||||
}
|
||||
|
||||
return outputPoints;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RiaPolyArcLineSampler::sampleLine( cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent )
|
||||
{
|
||||
cvf::Vec3d p1p2 = p2 - p1;
|
||||
|
||||
double p1p2Length = p1p2.length();
|
||||
|
||||
if ( m_isResamplingLines && p1p2Length > m_maxSamplingsInterval )
|
||||
{
|
||||
cvf::Vec3d tp1p2 = p1p2 / p1p2Length;
|
||||
double mdInc = m_maxSamplingsInterval;
|
||||
while ( mdInc < p1p2Length )
|
||||
{
|
||||
cvf::Vec3d ps = p1 + mdInc * tp1p2;
|
||||
m_points.push_back( ps );
|
||||
m_meshDs.push_back( m_totalMD + mdInc );
|
||||
mdInc += m_maxSamplingsInterval;
|
||||
}
|
||||
}
|
||||
|
||||
m_totalMD += p1p2Length;
|
||||
m_points.push_back( p2 );
|
||||
m_meshDs.push_back( m_totalMD );
|
||||
|
||||
( *endTangent ) = p1p2.getNormalized();
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RiaPolyArcLineSampler::sampleArc( cvf::Vec3d t1, cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent )
|
||||
{
|
||||
// Find arc CS
|
||||
RiaArcCurveCalculator CS_rad( p1, t1, p2 );
|
||||
|
||||
double radius = CS_rad.radius();
|
||||
cvf::Mat4d arcCS = CS_rad.arcCS();
|
||||
|
||||
double angleInc = m_maxSamplingsInterval / radius;
|
||||
|
||||
angleInc = angleInc < m_maxSamplingArcAngle ? angleInc : m_maxSamplingArcAngle; // Angle from 6 deg dogleg on 10 m
|
||||
|
||||
cvf::Vec3d C = CS_rad.center();
|
||||
cvf::Vec3d N = CS_rad.normal();
|
||||
|
||||
// Sample arc by
|
||||
// Rotate vector an increment, and transform to arc CS
|
||||
|
||||
double arcAngle = cvf::GeometryTools::getAngle( N, p1 - C, p2 - C );
|
||||
if ( arcAngle / angleInc > 5000 )
|
||||
{
|
||||
angleInc = arcAngle / 5000;
|
||||
}
|
||||
|
||||
for ( double angle = angleInc; angle < arcAngle; angle += angleInc )
|
||||
{
|
||||
cvf::Vec3d C_to_incP = cvf::Vec3d::X_AXIS;
|
||||
C_to_incP *= radius;
|
||||
C_to_incP.transformVector( cvf::Mat3d::fromRotation( cvf::Vec3d::Z_AXIS, angle ) );
|
||||
|
||||
C_to_incP.transformPoint( arcCS );
|
||||
|
||||
m_points.push_back( C_to_incP );
|
||||
m_meshDs.push_back( m_totalMD + angle * radius );
|
||||
}
|
||||
|
||||
m_totalMD += arcAngle * radius;
|
||||
m_points.push_back( p2 );
|
||||
m_meshDs.push_back( m_totalMD );
|
||||
|
||||
( *endTangent ) = CS_rad.endTangent();
|
||||
}
|
||||
@@ -0,0 +1,50 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2018- Equinor 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 "cvfVector3.h"
|
||||
|
||||
#include <vector>
|
||||
|
||||
class RiaPolyArcLineSampler
|
||||
{
|
||||
public:
|
||||
RiaPolyArcLineSampler( const cvf::Vec3d& startTangent, const std::vector<cvf::Vec3d>& lineArcEndPoints );
|
||||
|
||||
std::pair<std::vector<cvf::Vec3d>, std::vector<double>> sampledPointsAndMDs( double maxSampleInterval,
|
||||
bool isResamplingLines );
|
||||
|
||||
static std::vector<cvf::Vec3d> pointsWithoutDuplicates( const std::vector<cvf::Vec3d>& points );
|
||||
|
||||
private:
|
||||
void sampleLine( cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent );
|
||||
void sampleArc( cvf::Vec3d t1, cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent );
|
||||
void sampleSegment( cvf::Vec3d t1, cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent );
|
||||
|
||||
std::vector<cvf::Vec3d> m_points;
|
||||
std::vector<double> m_meshDs;
|
||||
|
||||
double m_maxSamplingsInterval;
|
||||
const double m_maxSamplingArcAngle = 0.07310818; // Angle from 6 deg dogleg on 10 m
|
||||
bool m_isResamplingLines;
|
||||
double m_totalMD;
|
||||
|
||||
cvf::Vec3d m_startTangent;
|
||||
std::vector<cvf::Vec3d> m_lineArcEndPoints;
|
||||
};
|
||||
@@ -0,0 +1,478 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2018- Equinor 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 "RiaSCurveCalculator.h"
|
||||
|
||||
#include "RiaOffshoreSphericalCoords.h"
|
||||
|
||||
#include "cvfMatrix4.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <iostream>
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
RiaSCurveCalculator::RiaSCurveCalculator( cvf::Vec3d p1,
|
||||
double azi1,
|
||||
double inc1,
|
||||
double rad1,
|
||||
cvf::Vec3d p2,
|
||||
double azi2,
|
||||
double inc2,
|
||||
double rad2 )
|
||||
: m_isCalculationOK( false )
|
||||
, m_p1( p1 )
|
||||
, m_p2( p2 )
|
||||
, m_firstArcEndpoint( p1 + 0.3 * ( p2 - p1 ) )
|
||||
, m_secondArcStartpoint( p1 + 0.6 * ( p2 - p1 ) )
|
||||
, m_r1( rad1 )
|
||||
, m_r2( rad2 )
|
||||
, m_ctrlPpointCurveStatus( NOT_SET )
|
||||
, m_solveStatus( NOT_SOLVED )
|
||||
{
|
||||
initializeByFinding_q1q2( p1, azi1, inc1, rad1, p2, azi2, inc2, rad2 );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
RiaSCurveCalculator::RiaSCurveCalculator( cvf::Vec3d p1, cvf::Vec3d q1, cvf::Vec3d p2, cvf::Vec3d q2 )
|
||||
: m_isCalculationOK( true )
|
||||
, m_p1( p1 )
|
||||
, m_p2( p2 )
|
||||
, m_ctrlPpointCurveStatus( NOT_SET )
|
||||
, m_solveStatus( NOT_SOLVED )
|
||||
{
|
||||
using Vec3d = cvf::Vec3d;
|
||||
bool isOk = true;
|
||||
m_isCalculationOK = true;
|
||||
|
||||
Vec3d tq1q2 = ( q2 - q1 ).getNormalized( &isOk ); // !ok means the control points are in the same place. Could
|
||||
// fallback to use only one circle segment + one line.
|
||||
m_isCalculationOK = m_isCalculationOK && isOk;
|
||||
Vec3d t1 = ( q1 - p1 ).getNormalized( &isOk ); // !ok means no tangent specified. Could fallback to use only one
|
||||
// circle segment + one line.
|
||||
m_isCalculationOK = m_isCalculationOK && isOk;
|
||||
Vec3d t2 = ( p2 - q2 ).getNormalized( &isOk ); // !ok means no tangent specified. Could fallback to use only one
|
||||
// circle segment + one line or only one straight line if both
|
||||
// tangents are missing
|
||||
m_isCalculationOK = m_isCalculationOK && isOk;
|
||||
|
||||
if ( !m_isCalculationOK )
|
||||
{
|
||||
m_ctrlPpointCurveStatus = FAILED_INPUT_OVERLAP;
|
||||
}
|
||||
|
||||
{
|
||||
Vec3d td1 = ( tq1q2 - t1 );
|
||||
double td1Length = td1.length();
|
||||
|
||||
if ( td1Length > 1e-10 )
|
||||
{
|
||||
td1 /= td1Length;
|
||||
m_c1 = q1 + ( ( q1 - p1 ).length() / ( td1 * ( -t1 ) ) ) * td1;
|
||||
m_r1 = ( m_c1 - p1 ).length();
|
||||
}
|
||||
else // both control points are along t1. First curve has infinite radius
|
||||
{
|
||||
m_c1 = cvf::Vec3d::UNDEFINED;
|
||||
m_r1 = std::numeric_limits<double>::infinity();
|
||||
|
||||
if ( m_ctrlPpointCurveStatus == NOT_SET )
|
||||
{
|
||||
m_ctrlPpointCurveStatus = OK_INFINITE_RADIUS1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
{
|
||||
Vec3d td2 = ( -tq1q2 + t2 );
|
||||
double td2Length = td2.length();
|
||||
|
||||
if ( td2Length > 1e-10 )
|
||||
{
|
||||
td2 /= td2Length;
|
||||
m_c2 = q2 + ( ( q2 - p2 ).length() / ( td2 * ( t2 ) ) ) * td2;
|
||||
m_r2 = ( m_c2 - p2 ).length();
|
||||
}
|
||||
else // both control points are along t2. Second curve has infinite radius
|
||||
{
|
||||
m_c2 = cvf::Vec3d::UNDEFINED;
|
||||
m_r2 = std::numeric_limits<double>::infinity();
|
||||
|
||||
if ( m_ctrlPpointCurveStatus == NOT_SET )
|
||||
{
|
||||
m_ctrlPpointCurveStatus = OK_INFINITE_RADIUS2;
|
||||
}
|
||||
else if ( m_ctrlPpointCurveStatus == OK_INFINITE_RADIUS1 )
|
||||
{
|
||||
m_ctrlPpointCurveStatus = OK_INFINITE_RADIUS12;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
m_firstArcEndpoint = q1 + ( q1 - p1 ).length() * tq1q2;
|
||||
m_secondArcStartpoint = q2 - ( q2 - p2 ).length() * tq1q2;
|
||||
|
||||
if ( ( ( q1 - p1 ).length() + ( q2 - p2 ).length() ) > ( q2 - q1 ).length() ) // first arc end and second arc start
|
||||
// is overlapping
|
||||
{
|
||||
m_ctrlPpointCurveStatus = FAILED_ARC_OVERLAP;
|
||||
m_isCalculationOK = false;
|
||||
}
|
||||
|
||||
if ( m_ctrlPpointCurveStatus == NOT_SET )
|
||||
{
|
||||
m_ctrlPpointCurveStatus = OK;
|
||||
}
|
||||
|
||||
// The Circle normals. Will be set to cvf::Vec3d::ZERO if undefined.
|
||||
|
||||
m_n1 = ( t1 ^ tq1q2 ).getNormalized();
|
||||
m_n2 = ( tq1q2 ^ t2 ).getNormalized();
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
RiaSCurveCalculator RiaSCurveCalculator::fromTangentsAndLength( cvf::Vec3d p1,
|
||||
double azi1,
|
||||
double inc1,
|
||||
double lengthToQ1,
|
||||
cvf::Vec3d p2,
|
||||
double azi2,
|
||||
double inc2,
|
||||
double lengthToQ2 )
|
||||
{
|
||||
cvf::Vec3d t1( RiaOffshoreSphericalCoords::unitVectorFromAziInc( azi1, inc1 ) );
|
||||
cvf::Vec3d t2( RiaOffshoreSphericalCoords::unitVectorFromAziInc( azi2, inc2 ) );
|
||||
|
||||
cvf::Vec3d Q1 = p1 + lengthToQ1 * t1;
|
||||
cvf::Vec3d Q2 = p2 - lengthToQ2 * t2;
|
||||
|
||||
RiaSCurveCalculator curveFromControlPoints( p1, Q1, p2, Q2 );
|
||||
|
||||
return curveFromControlPoints;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
/// Needs to calculate J^-1 * [R1_error, R2_error]
|
||||
/// | dR1_dq1 dR1_dq2 | 1 | dR2_dq2 -dR1_dq2 |
|
||||
/// J = | | J^-1 = ---------------------------------- | |
|
||||
/// | dR2_dq1 dR2_dq2 | dR1_dq1*dR2_dq2 - dR1_dq2*dR2_dq1 | -dR2_dq1 dR1_dq1 |
|
||||
///
|
||||
/// | q1_step | | R1_Error |
|
||||
/// | | = - J^-1 | |
|
||||
/// | q2_step | | R2_Error |
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void calculateNewStepsFromJacobi( double dR1_dq1,
|
||||
double dR1_dq2,
|
||||
double dR2_dq1,
|
||||
double dR2_dq2,
|
||||
double R1_error,
|
||||
double R2_error,
|
||||
double* newStepq1,
|
||||
double* newStepq2 )
|
||||
{
|
||||
double invJacobiScale = 1.0 / ( dR1_dq1 * dR2_dq2 - dR2_dq1 * dR1_dq2 );
|
||||
double invJacobi_R1q1 = invJacobiScale * dR2_dq2;
|
||||
double invJacobi_R1q2 = invJacobiScale * -dR1_dq2;
|
||||
double invJacobi_R2q1 = invJacobiScale * -dR2_dq1;
|
||||
double invJacobi_R2q2 = invJacobiScale * dR1_dq1;
|
||||
|
||||
( *newStepq1 ) = -( invJacobi_R1q1 * R1_error + invJacobi_R1q2 * R2_error );
|
||||
( *newStepq2 ) = -( invJacobi_R2q1 * R1_error + invJacobi_R2q2 * R2_error );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
bool isZeroCrossing( double newError, double oldError, double maxError )
|
||||
{
|
||||
if ( ( newError < -maxError && maxError < oldError ) || ( newError > maxError && -maxError > oldError ) )
|
||||
{
|
||||
return true;
|
||||
}
|
||||
|
||||
return false;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
/// Iterating with changing q1, q2 (lengths along tangent) to find solution with R1 = r1 and R2 = r2
|
||||
/// R1(q1, q2), R2(q1, q2)
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RiaSCurveCalculator::initializeByFinding_q1q2( cvf::Vec3d p1,
|
||||
double azi1,
|
||||
double inc1,
|
||||
double r1,
|
||||
cvf::Vec3d p2,
|
||||
double azi2,
|
||||
double inc2,
|
||||
double r2 )
|
||||
{
|
||||
// Algorithm options
|
||||
|
||||
const int maxIterations = 40;
|
||||
const double maxError = 0.01;
|
||||
const double maxStepSize = 1.0e9;
|
||||
const double maxLengthToQ = 1.0e10;
|
||||
bool enableBackstepping = false;
|
||||
//#define USE_JACOBI_UPDATE
|
||||
//#define DEBUG_OUTPUT_ON
|
||||
|
||||
// Needs the initial partial derivatives to see the direction of change
|
||||
// dR1/dq1, dR1/dq2, dR2/dq1, dR2/dq2
|
||||
// Selects a sensible point in the domain for the evaluation
|
||||
|
||||
double p1p2Length = ( p2 - p1 ).length();
|
||||
double delta = 0.01 * p1p2Length;
|
||||
double initialq1q2 = 0.2 * p1p2Length;
|
||||
double deltaPos = initialq1q2 + delta;
|
||||
|
||||
RiaSCurveCalculator ev_0 =
|
||||
RiaSCurveCalculator::fromTangentsAndLength( p1, azi1, inc1, initialq1q2, p2, azi2, inc2, initialq1q2 );
|
||||
|
||||
if ( ev_0.curveStatus() == RiaSCurveCalculator::OK_INFINITE_RADIUS12 )
|
||||
{
|
||||
*this = ev_0;
|
||||
this->m_solveStatus = CONVERGED;
|
||||
return;
|
||||
} // Todo: Handle infinite radius in one place
|
||||
|
||||
RiaSCurveCalculator ev_dq1 =
|
||||
RiaSCurveCalculator::fromTangentsAndLength( p1, azi1, inc1, deltaPos, p2, azi2, inc2, initialq1q2 );
|
||||
RiaSCurveCalculator ev_dq2 =
|
||||
RiaSCurveCalculator::fromTangentsAndLength( p1, azi1, inc1, initialq1q2, p2, azi2, inc2, deltaPos );
|
||||
|
||||
// Initial Jacobi
|
||||
double dR1_dq1 = ( ( r1 - ev_dq1.firstRadius() ) - ( r1 - ev_0.firstRadius() ) ) / delta;
|
||||
double dR2_dq2 = ( ( r2 - ev_dq2.secondRadius() ) - ( r2 - ev_0.secondRadius() ) ) / delta;
|
||||
|
||||
// Initial function value (error)
|
||||
double R1_error = r1 - ev_0.firstRadius();
|
||||
double R2_error = r2 - ev_0.secondRadius();
|
||||
|
||||
// First steps
|
||||
double q1Step = -R1_error / dR1_dq1;
|
||||
double q2Step = -R2_error / dR2_dq2;
|
||||
|
||||
#ifdef USE_JACOBI_UPDATE
|
||||
double dR1_dq2 = ( ( r1 - ev_dq2.firstRadius() ) - ( r1 - ev_0.firstRadius() ) ) / delta;
|
||||
double dR2_dq1 = ( ( r2 - ev_dq1.secondRadius() ) - ( r2 - ev_0.secondRadius() ) ) / delta;
|
||||
|
||||
calculateNewStepsFromJacobi( dR1_dq1, dR1_dq2, dR2_dq1, dR2_dq2, R1_error, R2_error, &q1Step, &q2Step );
|
||||
#endif
|
||||
|
||||
double q1 = initialq1q2;
|
||||
double q2 = initialq1q2;
|
||||
|
||||
#ifdef DEBUG_OUTPUT_ON
|
||||
std::cout << std::endl;
|
||||
std::cout << "Targets: R1, R2: " << r1 << " , " << r2 << std::endl;
|
||||
|
||||
std::cout << 0 << ": " << q1Step << " , " << q2Step << " : " << q1 << " , " << q2 << " | " << ev_0.isOk() << " : "
|
||||
<< ev_0.firstRadius() << " , " << ev_0.secondRadius() << " : " << R1_error << " , " << R2_error
|
||||
<< std::endl;
|
||||
#endif
|
||||
|
||||
SolveStatus solveResultStatus = NOT_SOLVED;
|
||||
|
||||
int backstepLevel = 0;
|
||||
int iteration = 1;
|
||||
for ( iteration = 1; iteration < maxIterations; ++iteration )
|
||||
{
|
||||
if ( fabs( q1Step ) > maxStepSize || fabs( q2Step ) > maxStepSize )
|
||||
{
|
||||
solveResultStatus = FAILED_MAX_TANGENT_STEP_REACHED;
|
||||
|
||||
break;
|
||||
}
|
||||
|
||||
std::string q1R1StepCorrMarker;
|
||||
std::string q2R2StepCorrMarker;
|
||||
|
||||
if ( q1 + q1Step < 0 )
|
||||
{
|
||||
q1Step = -0.9 * q1;
|
||||
q1R1StepCorrMarker = "*";
|
||||
}
|
||||
if ( q2 + q2Step < 0 )
|
||||
{
|
||||
q2Step = -0.9 * q2;
|
||||
q2R2StepCorrMarker = "*";
|
||||
}
|
||||
|
||||
q1 += q1Step;
|
||||
q2 += q2Step;
|
||||
|
||||
if ( fabs( q1 ) > maxLengthToQ || fabs( q2 ) > maxLengthToQ )
|
||||
{
|
||||
/// Max length along tangent reached
|
||||
solveResultStatus = FAILED_MAX_LENGTH_ALONG_TANGENT_REACHED;
|
||||
break;
|
||||
}
|
||||
|
||||
RiaSCurveCalculator ev_1 = RiaSCurveCalculator::fromTangentsAndLength( p1, azi1, inc1, q1, p2, azi2, inc2, q2 );
|
||||
|
||||
double R1_error_new = r1 - ev_1.firstRadius();
|
||||
double R2_error_new = r2 - ev_1.secondRadius();
|
||||
|
||||
#ifdef DEBUG_OUTPUT_ON
|
||||
std::cout << iteration << ": " << q1Step << q1R1StepCorrMarker << " , " << q2Step << q2R2StepCorrMarker << " : "
|
||||
<< q1 << " , " << q2 << " | " << ev_1.isOk() << " : " << ev_1.firstRadius() << " , "
|
||||
<< ev_1.secondRadius() << " : " << R1_error_new << " , " << R2_error_new;
|
||||
#endif
|
||||
|
||||
if ( ( fabs( R1_error_new ) < maxError || ev_1.curveStatus() == OK_INFINITE_RADIUS1 ) &&
|
||||
( fabs( R2_error_new ) < maxError || ev_1.curveStatus() == OK_INFINITE_RADIUS2 ) )
|
||||
{
|
||||
ev_0 = ev_1;
|
||||
|
||||
// Result ok !
|
||||
|
||||
solveResultStatus = CONVERGED;
|
||||
|
||||
#ifdef DEBUG_OUTPUT_ON
|
||||
std::cout << std::endl;
|
||||
#endif
|
||||
|
||||
break;
|
||||
}
|
||||
|
||||
if ( enableBackstepping ) // Experimental back-stepping
|
||||
{
|
||||
bool isZeroCrossingR1 = isZeroCrossing( R1_error_new, R1_error, maxError );
|
||||
bool isZeroCrossingR2 = isZeroCrossing( R2_error_new, R2_error, maxError );
|
||||
|
||||
if ( isZeroCrossingR2 || isZeroCrossingR1 )
|
||||
{
|
||||
q1 -= q1Step;
|
||||
q2 -= q2Step;
|
||||
|
||||
// if (isZeroCrossingR1)
|
||||
q1Step = 0.9 * q1Step * fabs( R1_error ) / ( fabs( R1_error_new ) + fabs( R1_error ) );
|
||||
// if (isZeroCrossingR2)
|
||||
q2Step = 0.9 * q2Step * fabs( R2_error ) / ( fabs( R2_error_new ) + fabs( R2_error ) );
|
||||
|
||||
++backstepLevel;
|
||||
|
||||
#ifdef DEBUG_OUTPUT_ON
|
||||
std::cout << " Backstep needed. " << std::endl;
|
||||
#endif
|
||||
|
||||
continue;
|
||||
}
|
||||
else
|
||||
{
|
||||
backstepLevel = 0;
|
||||
}
|
||||
}
|
||||
|
||||
#ifdef DEBUG_OUTPUT_ON
|
||||
std::cout << std::endl;
|
||||
#endif
|
||||
|
||||
#ifdef USE_JACOBI_UPDATE
|
||||
|
||||
/// Update Jacobi using Broyden
|
||||
// (R_error_n-Rerror_n-1) - Jn-1*dq
|
||||
// J_n = Jn-1 + --------------------------------- (dq)T
|
||||
// | dqn |^2
|
||||
//
|
||||
|
||||
double dR1_error = R1_error_new - R1_error;
|
||||
double dR2_error = R2_error_new - R2_error;
|
||||
R1_error = R1_error_new;
|
||||
R2_error = R2_error_new;
|
||||
|
||||
double stepNormScale = 1.0 / ( q1Step * q1Step + q2Step * q2Step );
|
||||
|
||||
dR1_dq1 = dR1_dq1 + stepNormScale * ( q1Step * ( dR1_error - q1Step * dR1_dq1 + q2Step * dR1_dq2 ) );
|
||||
dR1_dq2 = dR1_dq2 + stepNormScale * ( q2Step * ( dR1_error - q1Step * dR1_dq1 + q2Step * dR1_dq2 ) );
|
||||
dR2_dq1 = dR2_dq1 + stepNormScale * ( q1Step * ( dR2_error - q1Step * dR2_dq1 + q2Step * dR2_dq2 ) );
|
||||
dR2_dq2 = dR2_dq2 + stepNormScale * ( q2Step * ( dR2_error - q1Step * dR2_dq1 + q2Step * dR2_dq2 ) );
|
||||
|
||||
calculateNewStepsFromJacobi( dR1_dq1, dR1_dq2, dR2_dq1, dR2_dq2, R1_error, R2_error, &q1Step, &q2Step );
|
||||
|
||||
#else
|
||||
|
||||
dR1_dq1 = ( ( r1 - ev_1.firstRadius() ) - ( r1 - ev_0.firstRadius() ) ) / q1Step;
|
||||
dR2_dq2 = ( ( r2 - ev_1.secondRadius() ) - ( r2 - ev_0.secondRadius() ) ) / q2Step;
|
||||
|
||||
R1_error = R1_error_new;
|
||||
R2_error = R2_error_new;
|
||||
|
||||
q1Step = -R1_error / dR1_dq1;
|
||||
q2Step = -R2_error / dR2_dq2;
|
||||
|
||||
#endif
|
||||
|
||||
ev_0 = ev_1;
|
||||
}
|
||||
|
||||
*this = ev_0;
|
||||
if ( iteration >= maxIterations )
|
||||
{
|
||||
m_solveStatus = FAILED_MAX_ITERATIONS_REACHED;
|
||||
// Max iterations reached
|
||||
}
|
||||
else
|
||||
{
|
||||
m_solveStatus = solveResultStatus;
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RiaSCurveCalculator::dump() const
|
||||
{
|
||||
cvf::Vec3d v_C1 = firstCenter();
|
||||
cvf::Vec3d v_C2 = secondCenter();
|
||||
cvf::Vec3d v_N1 = firstNormal();
|
||||
cvf::Vec3d v_N2 = secondNormal();
|
||||
cvf::Vec3d v_P11 = firstArcEndpoint();
|
||||
cvf::Vec3d v_P22 = secondArcStartpoint();
|
||||
|
||||
std::cout << " P1: "
|
||||
<< "[ " << m_p1[0] << " " << m_p1[1] << " " << m_p1[2] << " " << std::endl;
|
||||
std::cout << " P11: "
|
||||
<< "[ " << v_P11[0] << " " << v_P11[1] << " " << v_P11[2] << " " << std::endl;
|
||||
std::cout << " P22: "
|
||||
<< "[ " << v_P22[0] << " " << v_P22[1] << " " << v_P22[2] << " " << std::endl;
|
||||
std::cout << " P2: "
|
||||
<< "[ " << m_p2[0] << " " << m_p2[1] << " " << m_p2[2] << " " << std::endl;
|
||||
std::cout << " C1: "
|
||||
<< "[ " << v_C1[0] << " " << v_C1[1] << " " << v_C1[2] << " " << std::endl;
|
||||
std::cout << " C2: "
|
||||
<< "[ " << v_C2[0] << " " << v_C2[1] << " " << v_C2[2] << " " << std::endl;
|
||||
std::cout << " N1: "
|
||||
<< "[ " << v_N1[0] << " " << v_N1[1] << " " << v_N1[2] << " " << std::endl;
|
||||
std::cout << " N2: "
|
||||
<< "[ " << v_N2[0] << " " << v_N2[1] << " " << v_N2[2] << " " << std::endl;
|
||||
std::cout << " R1: "
|
||||
<< "[ " << firstRadius() << " ]" << std::endl;
|
||||
std::cout << " R2: "
|
||||
<< "[ " << secondRadius() << " ]" << std::endl;
|
||||
std::cout << " CtrPointStatus: " << m_ctrlPpointCurveStatus << " SolveStatus: " << m_solveStatus << std::endl;
|
||||
}
|
||||
@@ -0,0 +1,98 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2018- Equinor 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 "cvfVector3.h"
|
||||
|
||||
class RiaSCurveCalculator
|
||||
{
|
||||
public:
|
||||
RiaSCurveCalculator( cvf::Vec3d p1, double azi1, double inc1, double r1, cvf::Vec3d p2, double azi2, double inc2, double r2 );
|
||||
|
||||
RiaSCurveCalculator( cvf::Vec3d p1, cvf::Vec3d q1, cvf::Vec3d p2, cvf::Vec3d q2 );
|
||||
|
||||
enum CurveStatus
|
||||
{
|
||||
NOT_SET,
|
||||
OK,
|
||||
OK_INFINITE_RADIUS1,
|
||||
OK_INFINITE_RADIUS2,
|
||||
OK_INFINITE_RADIUS12,
|
||||
FAILED_INPUT_OVERLAP,
|
||||
FAILED_ARC_OVERLAP
|
||||
};
|
||||
enum SolveStatus
|
||||
{
|
||||
NOT_SOLVED,
|
||||
CONVERGED,
|
||||
FAILED_MAX_LENGTH_ALONG_TANGENT_REACHED,
|
||||
FAILED_MAX_TANGENT_STEP_REACHED,
|
||||
FAILED_MAX_ITERATIONS_REACHED
|
||||
};
|
||||
|
||||
bool isOk() const { return m_isCalculationOK; }
|
||||
CurveStatus curveStatus() const { return m_ctrlPpointCurveStatus; }
|
||||
SolveStatus solveStatus() const { return m_solveStatus; }
|
||||
|
||||
cvf::Vec3d firstArcEndpoint() const { return m_firstArcEndpoint; }
|
||||
cvf::Vec3d secondArcStartpoint() const { return m_secondArcStartpoint; }
|
||||
cvf::Vec3d firstCenter() const { return m_c1; }
|
||||
cvf::Vec3d secondCenter() const { return m_c2; }
|
||||
cvf::Vec3d firstNormal() const { return m_n1; }
|
||||
cvf::Vec3d secondNormal() const { return m_n2; }
|
||||
double firstRadius() const { return m_r1; }
|
||||
double secondRadius() const { return m_r2; }
|
||||
|
||||
void dump() const;
|
||||
|
||||
static RiaSCurveCalculator fromTangentsAndLength( cvf::Vec3d p1,
|
||||
double azi1,
|
||||
double inc1,
|
||||
double lengthToQ1,
|
||||
cvf::Vec3d p2,
|
||||
double azi2,
|
||||
double inc2,
|
||||
double lengthToQ2 );
|
||||
|
||||
private:
|
||||
void initializeByFinding_q1q2( cvf::Vec3d p1,
|
||||
double azi1,
|
||||
double inc1,
|
||||
double r1,
|
||||
cvf::Vec3d p2,
|
||||
double azi2,
|
||||
double inc2,
|
||||
double r2 );
|
||||
|
||||
bool m_isCalculationOK;
|
||||
|
||||
CurveStatus m_ctrlPpointCurveStatus;
|
||||
SolveStatus m_solveStatus;
|
||||
|
||||
cvf::Vec3d m_p1;
|
||||
cvf::Vec3d m_p2;
|
||||
cvf::Vec3d m_firstArcEndpoint;
|
||||
cvf::Vec3d m_secondArcStartpoint;
|
||||
cvf::Vec3d m_c1;
|
||||
cvf::Vec3d m_c2;
|
||||
cvf::Vec3d m_n1;
|
||||
cvf::Vec3d m_n2;
|
||||
double m_r1;
|
||||
double m_r2;
|
||||
};
|
||||
@@ -0,0 +1,134 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2018- Equinor 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 "RiaWellPlanCalculator.h"
|
||||
#include "RiaArcCurveCalculator.h"
|
||||
#include "RiaOffshoreSphericalCoords.h"
|
||||
|
||||
#include "cvfGeometryTools.h"
|
||||
#include "cvfMatrix4.h"
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
RiaWellPlanCalculator::RiaWellPlanCalculator( const cvf::Vec3d& startTangent,
|
||||
const std::vector<cvf::Vec3d>& lineArcEndPoints )
|
||||
: m_startTangent( startTangent )
|
||||
, m_lineArcEndPoints( lineArcEndPoints )
|
||||
{
|
||||
if ( m_lineArcEndPoints.size() < 2 ) return;
|
||||
|
||||
WellPlanSegment segment = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
|
||||
|
||||
RiaOffshoreSphericalCoords startAziIncRad( m_startTangent );
|
||||
segment.inc = cvf::Math::toDegrees( startAziIncRad.inc() );
|
||||
segment.azi = cvf::Math::toDegrees( startAziIncRad.azi() );
|
||||
|
||||
segment.TVD = -lineArcEndPoints[0].z();
|
||||
segment.NS = lineArcEndPoints[0].y();
|
||||
segment.EW = lineArcEndPoints[0].x();
|
||||
|
||||
m_wpResult.push_back( segment );
|
||||
|
||||
cvf::Vec3d t2 = m_startTangent;
|
||||
|
||||
for ( size_t pIdx = 0; pIdx < m_lineArcEndPoints.size() - 1; ++pIdx )
|
||||
{
|
||||
addSegment( t2, m_lineArcEndPoints[pIdx], m_lineArcEndPoints[pIdx + 1], &t2 );
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RiaWellPlanCalculator::addSegment( cvf::Vec3d t1, cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent )
|
||||
{
|
||||
cvf::Vec3d p1p2 = p2 - p1;
|
||||
double xyLength = std::sqrt( p1p2.x() * p1p2.x() + p1p2.y() * p1p2.y() );
|
||||
double zLength = std::abs( p1p2.z() );
|
||||
|
||||
// We only show two decimals in the well plan anyway.
|
||||
if ( xyLength < 1.0e-2 && zLength < 1.0e-2 ) return;
|
||||
|
||||
if ( cvf::GeometryTools::getAngle( t1, p1p2 ) < 1e-5 )
|
||||
{
|
||||
addLineSegment( p1, p2, endTangent );
|
||||
}
|
||||
else // resample arc
|
||||
{
|
||||
addArcSegment( t1, p1, p2, endTangent );
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RiaWellPlanCalculator::addLineSegment( cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent )
|
||||
{
|
||||
WellPlanSegment segment = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
|
||||
|
||||
cvf::Vec3d p1p2 = p2 - p1;
|
||||
double length = p1p2.length();
|
||||
|
||||
segment.CL = length;
|
||||
segment.MD = m_wpResult.back().MD + length;
|
||||
|
||||
cvf::Vec3d tangent = p1p2 / length;
|
||||
|
||||
RiaOffshoreSphericalCoords aziIncRad( p1p2 );
|
||||
segment.inc = cvf::Math::toDegrees( aziIncRad.inc() );
|
||||
segment.azi = cvf::Math::toDegrees( aziIncRad.azi() );
|
||||
|
||||
segment.TVD = -p2.z();
|
||||
segment.NS = p2.y();
|
||||
segment.EW = p2.x();
|
||||
segment.dogleg = 0.0;
|
||||
segment.build = 0.0;
|
||||
segment.turn = 0.0;
|
||||
|
||||
m_wpResult.push_back( segment );
|
||||
*endTangent = tangent;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
void RiaWellPlanCalculator::addArcSegment( cvf::Vec3d t1, cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent )
|
||||
{
|
||||
WellPlanSegment segment = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
|
||||
|
||||
RiaArcCurveCalculator arcCalc( p1, t1, p2 );
|
||||
|
||||
segment.CL = arcCalc.arcLength();
|
||||
segment.MD = m_wpResult.back().MD + segment.CL;
|
||||
segment.inc = cvf::Math::toDegrees( arcCalc.endInclination() );
|
||||
segment.azi = cvf::Math::toDegrees( arcCalc.endAzimuth() );
|
||||
segment.TVD = -p2.z();
|
||||
segment.NS = p2.y();
|
||||
segment.EW = p2.x();
|
||||
segment.dogleg = cvf::Math::toDegrees( 30.0 / arcCalc.radius() );
|
||||
RiaOffshoreSphericalCoords startAziIncRad( t1 );
|
||||
double buildInRadsPrLength = ( arcCalc.endInclination() - startAziIncRad.inc() ) / arcCalc.arcLength();
|
||||
double turnInRadsPrLength = ( arcCalc.endAzimuth() - startAziIncRad.azi() ) / arcCalc.arcLength();
|
||||
segment.build = 30 * cvf::Math::toDegrees( buildInRadsPrLength );
|
||||
segment.turn = 30 * cvf::Math::toDegrees( turnInRadsPrLength );
|
||||
|
||||
m_wpResult.push_back( segment );
|
||||
|
||||
( *endTangent ) = arcCalc.endTangent();
|
||||
}
|
||||
@@ -0,0 +1,55 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2018- Equinor 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 "cvfVector3.h"
|
||||
|
||||
#include <vector>
|
||||
|
||||
class RiaWellPlanCalculator
|
||||
{
|
||||
public:
|
||||
RiaWellPlanCalculator( const cvf::Vec3d& startTangent, const std::vector<cvf::Vec3d>& lineArcEndPoints );
|
||||
|
||||
struct WellPlanSegment
|
||||
{
|
||||
double MD;
|
||||
double CL;
|
||||
double inc;
|
||||
double azi;
|
||||
double TVD;
|
||||
double NS;
|
||||
double EW;
|
||||
double dogleg;
|
||||
double build;
|
||||
double turn;
|
||||
};
|
||||
|
||||
const std::vector<WellPlanSegment>& wellPlan() const { return m_wpResult; }
|
||||
|
||||
private:
|
||||
void addSegment( cvf::Vec3d t1, cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent );
|
||||
void addLineSegment( cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent );
|
||||
void addArcSegment( cvf::Vec3d t1, cvf::Vec3d p1, cvf::Vec3d p2, cvf::Vec3d* endTangent );
|
||||
|
||||
cvf::Vec3d m_startTangent;
|
||||
std::vector<cvf::Vec3d> m_lineArcEndPoints;
|
||||
|
||||
std::vector<WellPlanSegment> m_wpResult;
|
||||
};
|
||||
Reference in New Issue
Block a user