#6544 Improve stress anisotropy time lapse calculation.

This commit is contained in:
Kristian Bendiksen
2020-09-28 15:17:06 +02:00
parent 0ee12f5166
commit 18c6d648e6
3 changed files with 116 additions and 0 deletions

View File

@@ -130,3 +130,112 @@ RigFemScalarResultFrames* RigFemPartResultCalculatorStressAnisotropy::calculate(
RigFemScalarResultFrames* requestedStress = m_resultCollection->findOrLoadScalarResult( partIndex, resVarAddr );
return requestedStress;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigFemScalarResultFrames*
RigFemPartResultCalculatorStressAnisotropy::calculateTimeLapse( int partIndex, const RigFemResultAddress& resVarAddr )
{
CVF_ASSERT( isMatching( resVarAddr ) );
caf::ProgressInfo frameCountProgress( m_resultCollection->frameCount() * 4, "" );
frameCountProgress.setProgressDescription(
"Calculating " + QString::fromStdString( resVarAddr.fieldName + ": " + resVarAddr.componentName ) );
frameCountProgress.setNextProgressIncrement( m_resultCollection->frameCount() );
RigFemScalarResultFrames* s1Frames =
m_resultCollection->findOrLoadScalarResult( partIndex,
RigFemResultAddress( resVarAddr.resultPosType,
resVarAddr.fieldName,
"S1" ) );
frameCountProgress.incrementProgress();
frameCountProgress.setNextProgressIncrement( m_resultCollection->frameCount() );
RigFemScalarResultFrames* s2Frames =
m_resultCollection->findOrLoadScalarResult( partIndex,
RigFemResultAddress( resVarAddr.resultPosType,
resVarAddr.fieldName,
"S2" ) );
frameCountProgress.incrementProgress();
frameCountProgress.setNextProgressIncrement( m_resultCollection->frameCount() );
RigFemScalarResultFrames* s3Frames =
m_resultCollection->findOrLoadScalarResult( partIndex,
RigFemResultAddress( resVarAddr.resultPosType,
resVarAddr.fieldName,
"S3" ) );
int baseTimeStep = resVarAddr.timeLapseBaseFrameIdx;
RigFemScalarResultFrames* s12Frames =
m_resultCollection->createScalarResult( partIndex,
RigFemResultAddress( resVarAddr.resultPosType,
resVarAddr.fieldName,
"SA12",
resVarAddr.timeLapseBaseFrameIdx ) );
RigFemScalarResultFrames* s13Frames =
m_resultCollection->createScalarResult( partIndex,
RigFemResultAddress( resVarAddr.resultPosType,
resVarAddr.fieldName,
"SA13",
resVarAddr.timeLapseBaseFrameIdx ) );
RigFemScalarResultFrames* s23Frames =
m_resultCollection->createScalarResult( partIndex,
RigFemResultAddress( resVarAddr.resultPosType,
resVarAddr.fieldName,
"SA23",
resVarAddr.timeLapseBaseFrameIdx ) );
frameCountProgress.incrementProgress();
frameCountProgress.setNextProgressIncrement( 1u );
float inf = std::numeric_limits<float>::infinity();
int frameCount = s1Frames->frameCount();
for ( int fIdx = 0; fIdx < frameCount; ++fIdx )
{
const std::vector<float>& s1t = s1Frames->frameData( fIdx );
const std::vector<float>& s2t = s2Frames->frameData( fIdx );
const std::vector<float>& s3t = s3Frames->frameData( fIdx );
const std::vector<float>& s1b = s1Frames->frameData( baseTimeStep );
const std::vector<float>& s2b = s2Frames->frameData( baseTimeStep );
const std::vector<float>& s3b = s3Frames->frameData( baseTimeStep );
std::vector<float>& s12 = s12Frames->frameData( fIdx );
std::vector<float>& s13 = s13Frames->frameData( fIdx );
std::vector<float>& s23 = s23Frames->frameData( fIdx );
size_t valCount = s1t.size();
s12.resize( valCount, 0.0 );
s13.resize( valCount, 0.0 );
s23.resize( valCount, 0.0 );
#pragma omp parallel for schedule( dynamic )
for ( long vIdx = 0; vIdx < static_cast<long>( valCount ); ++vIdx )
{
if ( fIdx != baseTimeStep )
{
double diffS1 = s1t[vIdx] - s1b[vIdx];
double diffS2 = s2t[vIdx] - s2b[vIdx];
double diffS3 = s3t[vIdx] - s3b[vIdx];
if ( diffS1 + diffS2 != 0.0 )
s12[vIdx] = 2.0 * ( diffS1 - diffS2 ) / ( diffS1 + diffS2 );
else
s12[vIdx] = inf;
if ( diffS1 + diffS3 != 0.0 )
s13[vIdx] = 2.0 * ( diffS1 - diffS3 ) / ( diffS1 + diffS3 );
else
s13[vIdx] = inf;
if ( diffS2 + diffS3 != 0.0 )
s23[vIdx] = 2.0 * ( diffS2 - diffS3 ) / ( diffS2 + diffS3 );
else
s23[vIdx] = inf;
}
}
frameCountProgress.incrementProgress();
}
RigFemScalarResultFrames* requestedStress = m_resultCollection->findOrLoadScalarResult( partIndex, resVarAddr );
return requestedStress;
}

View File

@@ -34,4 +34,5 @@ public:
virtual ~RigFemPartResultCalculatorStressAnisotropy();
bool isMatching( const RigFemResultAddress& resVarAddr ) const override;
RigFemScalarResultFrames* calculate( int partIndex, const RigFemResultAddress& resVarAddr ) override;
RigFemScalarResultFrames* calculateTimeLapse( int partIndex, const RigFemResultAddress& resVarAddr );
};

View File

@@ -22,6 +22,7 @@
#include "RigFemPartCollection.h"
#include "RigFemPartResultCalculatorGamma.h"
#include "RigFemPartResultCalculatorNormalized.h"
#include "RigFemPartResultCalculatorStressAnisotropy.h"
#include "RigFemPartResultsCollection.h"
#include "RigFemResultAddress.h"
#include "RigFemScalarResultFrames.h"
@@ -61,10 +62,15 @@ RigFemScalarResultFrames* RigFemPartResultCalculatorTimeLapse::calculate( int
{
CVF_ASSERT( resVarAddr.isTimeLapse() );
RigFemPartResultCalculatorStressAnisotropy anisotropyCalculator( *m_resultCollection );
if ( resVarAddr.fieldName == "Gamma" )
{
return calculateGammaTimeLapse( partIndex, resVarAddr );
}
else if ( anisotropyCalculator.isMatching( resVarAddr ) )
{
return anisotropyCalculator.calculateTimeLapse( partIndex, resVarAddr );
}
else
{
return calculateTimeLapse( partIndex, resVarAddr );