OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

Displacement curves fault faces (#10243)

Browse Source
This commit is contained in:
Jørgen Herje
2023-05-22 11:37:57 +02:00
committed by GitHub
parent 9a219ddcb2
commit b83fe73395
33 changed files with 1235 additions and 47 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
ApplicationLibCode/GeoMech/GeoMechDataModel/CMakeLists.txt
View File

@@ -104,6 +104,8 @@ add_library(
RigFemPartResultCalculatorKIndices.cpp
RimGeoMechGeometrySelectionItem.h
RimGeoMechGeometrySelectionItem.cpp
RigFemPartResultCalculatorNodalDisplacement.h
RigFemPartResultCalculatorNodalDisplacement.cpp
)
target_include_directories(${PROJECT_NAME} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})

6
ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPart.cpp
View File

@@ -477,16 +477,16 @@ cvf::BoundingBox RigFemPart::boundingBox() const
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigFemPart::findIntersectingCells( const cvf::BoundingBox& inputBB, std::vector<size_t>* elementIndices ) const
void RigFemPart::findIntersectingElementIndices( const cvf::BoundingBox& inputBB, std::vector<size_t>* elementIndices ) const
{
ensureIntersectionSearchTreeIsBuilt();
findIntersectingCellsWithExistingSearchTree( inputBB, elementIndices );
findIntersectingElementsWithExistingSearchTree( inputBB, elementIndices );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigFemPart::findIntersectingCellsWithExistingSearchTree( const cvf::BoundingBox& inputBB, std::vector<size_t>* elementIndices ) const
void RigFemPart::findIntersectingElementsWithExistingSearchTree( const cvf::BoundingBox& inputBB, std::vector<size_t>* elementIndices ) const
{
CVF_ASSERT( m_elementSearchTree.notNull() );
m_elementSearchTree->findIntersections( inputBB, elementIndices );

4
ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPart.h
View File

@@ -83,8 +83,8 @@ public:
cvf::BoundingBox boundingBox() const;
float characteristicElementSize() const;
const std::vector<int>& possibleGridCornerElements() const;
void findIntersectingCells( const cvf::BoundingBox& inputBB, std::vector<size_t>* elementIndices ) const;
void findIntersectingCellsWithExistingSearchTree( const cvf::BoundingBox& inputBB, std::vector<size_t>* elementIndices ) const;
void findIntersectingElementIndices( const cvf::BoundingBox& inputBB, std::vector<size_t>* elementIndices ) const;
void findIntersectingElementsWithExistingSearchTree( const cvf::BoundingBox& inputBB, std::vector<size_t>* elementIndices ) const;
void ensureIntersectionSearchTreeIsBuilt() const;

26
ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartCollection.cpp
View File

@@ -127,7 +127,7 @@ cvf::BoundingBox RigFemPartCollection::boundingBox() const
}
//--------------------------------------------------------------------------------------------------
/// convert from global element index to part and part-local index
/// Get part id and local part element index from global element index
//--------------------------------------------------------------------------------------------------
std::pair<int, size_t> RigFemPartCollection::partIdAndElementIndex( size_t globalIndex ) const
{
@@ -147,7 +147,7 @@ std::pair<int, size_t> RigFemPartCollection::partIdAndElementIndex( size_t globa
}
//--------------------------------------------------------------------------------------------------
/// convert from global element index to part and part-local index
/// Get part and local part element index from global element index
//--------------------------------------------------------------------------------------------------
std::pair<const RigFemPart*, size_t> RigFemPartCollection::partAndElementIndex( size_t globalIndex ) const
{
@@ -156,13 +156,31 @@ std::pair<const RigFemPart*, size_t> RigFemPartCollection::partAndElementIndex(
}
//--------------------------------------------------------------------------------------------------
/// convert from part and part-local index to global index
/// Get global element index from part id and local part element index
//--------------------------------------------------------------------------------------------------
size_t RigFemPartCollection::globalIndex( int partId, size_t localIndex ) const
size_t RigFemPartCollection::globalElementIndex( int partId, size_t localIndex ) const
{
return localIndex + m_partElementOffset[partId];
}
//--------------------------------------------------------------------------------------------------
/// Find intersecting global element indexes for a given bounding box
//--------------------------------------------------------------------------------------------------
void RigFemPartCollection::findIntersectingGlobalElementIndices( const cvf::BoundingBox& intersectingBB,
std::vector<size_t>* intersectedGlobalElementIndices ) const
{
for ( const auto& part : m_femParts )
{
std::vector<size_t> foundElements;
part->findIntersectingElementIndices( intersectingBB, &foundElements );
for ( const auto& foundElement : foundElements )
{
const size_t globalIdx = globalElementIndex( part->elementPartId(), foundElement );
intersectedGlobalElementIndices->push_back( globalIdx );
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------

5
ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartCollection.h
View File

@@ -40,9 +40,10 @@ public:
std::pair<int, size_t> partIdAndElementIndex( size_t globalIndex ) const;
std::pair<const RigFemPart*, size_t> partAndElementIndex( size_t globalIndex ) const;
size_t globalIndex( int partId, size_t localIndex ) const;
size_t globalElementIndex( int partId, size_t localIndex ) const;
std::pair<int, size_t> partIdAndNodeIndex( size_t globalNodeIndex ) const;
void findIntersectingGlobalElementIndices( const cvf::BoundingBox& intersectingBB,
std::vector<size_t>* intersectedGlobalElementIndices ) const;
int nodeIdxFromElementNodeResultIdx( size_t globalResultIdx ) const;

2
ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorCompaction.cpp
View File

@@ -164,7 +164,7 @@ void findReferenceElementForNode( const RigFemPart& part, size_t nodeIdx, size_t
bb.add( p2 );
std::vector<size_t> refElementCandidates;
part.findIntersectingCells( bb, &refElementCandidates );
part.findIntersectingElementIndices( bb, &refElementCandidates );
const RigFemPartGrid* grid = part.getOrCreateStructGrid();

128
ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorNodalDisplacement.cpp Normal file
View File

@@ -0,0 +1,128 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2023- 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 "RigFemPartResultCalculatorNodalDisplacement.h"
#include "RigFemPart.h"
#include "RigFemPartCollection.h"
#include "RigFemPartResultsCollection.h"
#include "RigFemResultAddress.h"
#include "RigFemScalarResultFrames.h"
#include "RigHexGradientTools.h"
#include "cafProgressInfo.h"
#include <QString>
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigFemPartResultCalculatorNodalDisplacement::RigFemPartResultCalculatorNodalDisplacement( RigFemPartResultsCollection& collection )
: RigFemPartResultCalculator( collection )
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigFemPartResultCalculatorNodalDisplacement::~RigFemPartResultCalculatorNodalDisplacement()
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartResultCalculatorNodalDisplacement::isMatching( const RigFemResultAddress& resVarAddr ) const
{
return ( ( resVarAddr.resultPosType == RIG_NODAL ) && ( resVarAddr.fieldName == "U" ) && ( resVarAddr.componentName == "U_LENGTH" ) );
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigFemScalarResultFrames* RigFemPartResultCalculatorNodalDisplacement::calculate( int partIndex, const RigFemResultAddress& resVarAddr )
{
CVF_ASSERT( resVarAddr.fieldName == "U" );
CVF_ASSERT( resVarAddr.componentName == "U_LENGTH" );
// Why calculating by 5?
caf::ProgressInfo stepCountProgress( m_resultCollection->timeStepCount() * 5, "" );
stepCountProgress.setProgressDescription( "Calculating vector length: " +
QString::fromStdString( resVarAddr.fieldName + ": " + resVarAddr.componentName ) );
stepCountProgress.setNextProgressIncrement( m_resultCollection->timeStepCount() );
const RigFemScalarResultFrames* srcFramesU1 =
m_resultCollection->findOrLoadScalarResult( partIndex, RigFemResultAddress( RIG_NODAL, resVarAddr.fieldName, "U1" ) );
stepCountProgress.incrementProgress();
stepCountProgress.setNextProgressIncrement( m_resultCollection->timeStepCount() );
const RigFemScalarResultFrames* srcFramesU2 =
m_resultCollection->findOrLoadScalarResult( partIndex, RigFemResultAddress( RIG_NODAL, resVarAddr.fieldName, "U2" ) );
stepCountProgress.incrementProgress();
stepCountProgress.setNextProgressIncrement( m_resultCollection->timeStepCount() );
const RigFemScalarResultFrames* srcFramesU3 =
m_resultCollection->findOrLoadScalarResult( partIndex, RigFemResultAddress( RIG_NODAL, resVarAddr.fieldName, "U3" ) );
stepCountProgress.incrementProgress();
stepCountProgress.setNextProgressIncrement( m_resultCollection->timeStepCount() );
CVF_ASSERT( srcFramesU1->timeStepCount() == srcFramesU2->timeStepCount() && srcFramesU2->timeStepCount() == srcFramesU3->timeStepCount() );
RigFemResultAddress dstVarAddrULength( RIG_NODAL, resVarAddr.fieldName, resVarAddr.componentName );
RigFemScalarResultFrames* dstFramesULength = m_resultCollection->createScalarResult( partIndex, dstVarAddrULength );
stepCountProgress.incrementProgress();
const RigFemPart* femPart = m_resultCollection->parts()->part( partIndex );
constexpr float inf = std::numeric_limits<float>::infinity();
const int timeStepCount = m_resultCollection->timeStepCount();
for ( int tIdx = 0; tIdx < timeStepCount; ++tIdx )
{
const int frameCount = m_resultCollection->frameCount( tIdx );
for ( int fIdx = 0; fIdx < frameCount; ++fIdx )
{
std::vector<float>& dstFrameDataULength = dstFramesULength->frameData( tIdx, fIdx );
const std::vector<float>& srcFrameDataU1 = srcFramesU1->frameData( tIdx, fIdx );
const std::vector<float>& srcFrameDataU2 = srcFramesU2->frameData( tIdx, fIdx );
const std::vector<float>& srcFrameDataU3 = srcFramesU3->frameData( tIdx, fIdx );
// Create empty results if:
// - One or more vector components are missing
// - Vector components are not of equal size
if ( srcFrameDataU1.empty() || srcFrameDataU2.empty() || srcFrameDataU3.empty() ||
( srcFrameDataU1.size() != srcFrameDataU2.size() || srcFrameDataU2.size() != srcFrameDataU3.size() ) )
{
dstFrameDataULength.clear();
continue;
}
dstFrameDataULength.resize( srcFrameDataU1.size(), 0.0 );
#pragma omp parallel for
for ( int i = 0; i < static_cast<int>( dstFrameDataULength.size() ); ++i )
{
const float u1 = srcFrameDataU1[i];
const float u2 = srcFrameDataU2[i];
const float u3 = srcFrameDataU3[i];
const float uLength = sqrt( u1 * u1 + u2 * u2 + u3 * u3 );
dstFrameDataULength[i] = uLength;
}
}
stepCountProgress.incrementProgress();
}
return dstFramesULength;
}

37
ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultCalculatorNodalDisplacement.h Normal file
View File

@@ -0,0 +1,37 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2023- 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 "RigFemPartResultCalculator.h"
class RigFemPartResultsCollection;
class RigFemScalarResultFrames;
class RigFemResultAddress;
//==================================================================================================
///
//==================================================================================================
class RigFemPartResultCalculatorNodalDisplacement : public RigFemPartResultCalculator
{
public:
explicit RigFemPartResultCalculatorNodalDisplacement( RigFemPartResultsCollection& collection );
~RigFemPartResultCalculatorNodalDisplacement() override;
bool isMatching( const RigFemResultAddress& resVarAddr ) const override;
RigFemScalarResultFrames* calculate( int partIndex, const RigFemResultAddress& resVarAddr ) override;
};

22
ApplicationLibCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp
View File

@@ -41,6 +41,7 @@
#include "RigFemPartResultCalculatorKIndices.h"
#include "RigFemPartResultCalculatorMudWeightWindow.h"
#include "RigFemPartResultCalculatorNE.h"
#include "RigFemPartResultCalculatorNodalDisplacement.h"
#include "RigFemPartResultCalculatorNodalGradients.h"
#include "RigFemPartResultCalculatorNormalSE.h"
#include "RigFemPartResultCalculatorNormalST.h"
@@ -197,6 +198,7 @@ RigFemPartResultsCollection::RigFemPartResultsCollection( RifGeoMechReaderInterf
m_resultCalculators.push_back( std::unique_ptr<RigFemPartResultCalculator>( new RigFemPartResultCalculatorShearSlipIndicator( *this ) ) );
m_resultCalculators.push_back( std::unique_ptr<RigFemPartResultCalculator>( new RigFemPartResultCalculatorFormationIndices( *this ) ) );
m_resultCalculators.push_back( std::unique_ptr<RigFemPartResultCalculator>( new RigFemPartResultCalculatorKIndices( *this ) ) );
m_resultCalculators.push_back( std::unique_ptr<RigFemPartResultCalculator>( new RigFemPartResultCalculatorNodalDisplacement( *this ) ) );
}
//--------------------------------------------------------------------------------------------------
@@ -452,16 +454,7 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::findOrLoadScalarResult( i
currentFrames->enableAsSingleStepResult();
currentFrames->frameData( 0, 0 ).swap( values );
}
frames = m_femPartResults[partIndex]->findScalarResult( resVarAddr );
if ( frames )
{
return frames;
}
else
{
return m_femPartResults[partIndex]->createScalarResult( resVarAddr );
}
return m_femPartResults[partIndex]->createScalarResult( resVarAddr );
}
// We need to read the data as bulk fields, and populate the correct scalar caches
@@ -567,6 +560,7 @@ std::map<std::string, std::vector<std::string>> RigFemPartResultsCollection::sca
if ( resPos == RIG_NODAL )
{
fieldCompNames = m_readerInterface->scalarNodeFieldAndComponentNames();
if ( fieldCompNames.contains( "U" ) ) fieldCompNames["U"].push_back( "U_LENGTH" );
fieldCompNames["POR-Bar"];
fieldCompNames[FIELD_NAME_COMPACTION];
}
@@ -650,14 +644,14 @@ std::map<std::string, std::vector<std::string>> RigFemPartResultsCollection::sca
fieldCompNames["MUD-WEIGHT"].push_back( "UMWL" );
fieldCompNames["MUD-WEIGHT"].push_back( "LMWL" );
if ( fieldCompNames.count( "LE" ) > 0 )
if ( fieldCompNames.contains( "LE" ) )
{
fieldCompNames["LE"].push_back( "LE1" );
fieldCompNames["LE"].push_back( "LE2" );
fieldCompNames["LE"].push_back( "LE3" );
}
if ( fieldCompNames.count( "PE" ) > 0 )
if ( fieldCompNames.contains( "PE" ) )
{
fieldCompNames["PE"].push_back( "PE1" );
fieldCompNames["PE"].push_back( "PE2" );
@@ -756,14 +750,14 @@ std::map<std::string, std::vector<std::string>> RigFemPartResultsCollection::sca
fieldCompNames["MUD-WEIGHT"].push_back( "UMWL" );
fieldCompNames["MUD-WEIGHT"].push_back( "LMWL" );
if ( fieldCompNames.count( "LE" ) > 0 )
if ( fieldCompNames.contains( "LE" ) )
{
fieldCompNames["LE"].push_back( "LE1" );
fieldCompNames["LE"].push_back( "LE2" );
fieldCompNames["LE"].push_back( "LE3" );
}
if ( fieldCompNames.count( "PE" ) > 0 )
if ( fieldCompNames.contains( "PE" ) )
{
fieldCompNames["PE"].push_back( "PE1" );
fieldCompNames["PE"].push_back( "PE2" );