diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp index a0d003f235..47bb7dfb8a 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp @@ -101,8 +101,9 @@ RigFemPartResultsCollection::RigFemPartResultsCollection( RifGeoMechReaderInterf femPartResult->initResultSteps( filteredStepNames ); } - m_cohesion = 10.0; - m_frictionAngleRad = cvf::Math::toRadians( 30.0 ); + m_cohesion = 10.0; + m_frictionAngleRad = cvf::Math::toRadians( 30.0 ); + m_normalizationAirGap = 0.0; } //-------------------------------------------------------------------------------------------------- @@ -735,6 +736,10 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateTimeLapseResult( frameCountProgress.incrementProgress(); calculateGammaFromFrames( partIndex, srcDataFrames, srcPORDataFrames, dstDataFrames, &frameCountProgress ); + if ( resVarAddr.normalizeByHydrostaticPressure() ) + { + dstDataFrames = calculateNormalizedResult( partIndex, resVarAddr ); + } return dstDataFrames; } @@ -1203,6 +1208,77 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateNodalGradients( return requestedGradient; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemScalarResultFrames* RigFemPartResultsCollection::calculateNormalizedResult( int partIndex, + const RigFemResultAddress& resVarAddr ) +{ + CVF_ASSERT( resVarAddr.normalizeByHydrostaticPressure() && isNormalizableResult( resVarAddr ) ); + + RigFemResultAddress unscaledResult = resVarAddr; + unscaledResult.normalizedByHydrostaticPressure = false; + + caf::ProgressInfo frameCountProgress( this->frameCount() * 3, "Calculating Normalized Result" ); + + RigFemScalarResultFrames* srcDataFrames = nullptr; + RigFemScalarResultFrames* dstDataFrames = nullptr; + + { + auto task = frameCountProgress.task( "Loading Unscaled Result", this->frameCount() ); + srcDataFrames = this->findOrLoadScalarResult( partIndex, unscaledResult ); + if ( !srcDataFrames ) return nullptr; + } + { + auto task = frameCountProgress.task( "Creating Space for Normalized Result", this->frameCount() ); + dstDataFrames = m_femPartResults[partIndex]->createScalarResult( RigFemResultAddress( resVarAddr ) ); + if ( !dstDataFrames ) return nullptr; + } + + frameCountProgress.setNextProgressIncrement( 1u ); + + const RigFemPart* femPart = m_femParts->part( partIndex ); + float inf = std::numeric_limits::infinity(); + int elmNodeCount = femPart->elementCount(); + const std::vector& nodeCoords = femPart->nodes().coordinates; + + int frameCount = srcDataFrames->frameCount(); + for ( int fIdx = 0; fIdx < frameCount; ++fIdx ) + { + const std::vector& srcFrameData = srcDataFrames->frameData( fIdx ); + std::vector& dstFrameData = dstDataFrames->frameData( fIdx ); + + size_t resultCount = srcFrameData.size(); + dstFrameData.resize( resultCount ); + + if ( resVarAddr.resultPosType == RIG_ELEMENT_NODAL ) + { +#pragma omp parallel for schedule( dynamic ) + for ( int elmNodeResIdx = 0; elmNodeResIdx < resultCount; ++elmNodeResIdx ) + { + const int nodeIdx = femPart->nodeIdxFromElementNodeResultIdx( elmNodeResIdx ); + const cvf::Vec3f node = nodeCoords[nodeIdx]; + double tvdRKB = std::abs( node.z() ) + m_normalizationAirGap; + double hydrostaticPressure = RiaWellLogUnitTools::hydrostaticPorePressureBar( tvdRKB ); + dstFrameData[elmNodeResIdx] = srcFrameData[elmNodeResIdx] / hydrostaticPressure; + } + } + else if ( resVarAddr.resultPosType == RIG_NODAL ) + { +#pragma omp parallel for schedule( dynamic ) + for ( int nodeResIdx = 0; nodeResIdx < resultCount; ++nodeResIdx ) + { + const cvf::Vec3f node = nodeCoords[nodeResIdx]; + double tvdRKB = std::abs( node.z() ) + m_normalizationAirGap; + double hydrostaticPressure = RiaWellLogUnitTools::hydrostaticPorePressureBar( tvdRKB ); + dstFrameData[nodeResIdx] = srcFrameData[nodeResIdx] / hydrostaticPressure; + } + } + frameCountProgress.incrementProgress(); + } + return dstDataFrames; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -2328,6 +2404,11 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateDerivedResult( i return calculateTimeLapseResult( partIndex, resVarAddr ); } + if ( resVarAddr.normalizeByHydrostaticPressure() ) + { + return calculateNormalizedResult( partIndex, resVarAddr ); + } + if ( resVarAddr.resultPosType == RIG_ELEMENT_NODAL_FACE ) { if ( resVarAddr.componentName == "Pazi" || resVarAddr.componentName == "Pinc" ) @@ -2941,6 +3022,69 @@ std::vector return addresses; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::set RigFemPartResultsCollection::normalizedResults() +{ + std::set validFields = {"SE", "ST"}; + std::set validComponents = {"S11", "S22", "S33", "S12", "S13", "S23", "S1", "S2", "S3", "Q"}; + + std::set results; + for ( auto field : validFields ) + { + for ( auto component : validComponents ) + { + results.insert( RigFemResultAddress( RIG_ELEMENT_NODAL, + field, + component, + RigFemResultAddress::allTimeLapsesValue(), + -1, + true ) ); + } + } + results.insert( + RigFemResultAddress( RIG_ELEMENT_NODAL, "SE", "SEM", RigFemResultAddress::allTimeLapsesValue(), -1, true ) ); + results.insert( + RigFemResultAddress( RIG_ELEMENT_NODAL, "ST", "STM", RigFemResultAddress::allTimeLapsesValue(), -1, true ) ); + results.insert( RigFemResultAddress( RIG_NODAL, "POR-Bar", "", RigFemResultAddress::allTimeLapsesValue(), -1, true ) ); + results.insert( + RigFemResultAddress( RIG_ELEMENT_NODAL, "POR-Bar", "", RigFemResultAddress::allTimeLapsesValue(), -1, true ) ); + + return results; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +bool RigFemPartResultsCollection::isNormalizableResult( const RigFemResultAddress& result ) +{ + for ( auto normRes : normalizedResults() ) + { + if ( normRes.resultPosType == result.resultPosType && normRes.fieldName == result.fieldName && + normRes.componentName == result.componentName ) + { + return true; + } + } + return false; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RigFemPartResultsCollection::setNormalizationAirGap( double normalizationAirGap ) +{ + if ( std::abs( m_normalizationAirGap - normalizationAirGap ) > 1.0e-8 ) + { + for ( auto result : normalizedResults() ) + { + this->deleteResult( result ); + } + } + m_normalizationAirGap = normalizationAirGap; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h index 7f56ba80aa..b12a7f2ea2 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h @@ -121,6 +121,10 @@ public: static std::vector tensorComponentAddresses( const RigFemResultAddress& resVarAddr ); static std::vector tensorPrincipalComponentAdresses( const RigFemResultAddress& resVarAddr ); + static std::set normalizedResults(); + static bool isNormalizableResult( const RigFemResultAddress& result ); + + void setNormalizationAirGap( double normalizationAirGap ); private: RigFemScalarResultFrames* findOrLoadScalarResult( int partIndex, const RigFemResultAddress& resVarAddr ); @@ -159,6 +163,7 @@ private: RigFemScalarResultFrames* calculateFormationIndices( int partIndex, const RigFemResultAddress& resVarAddr ); RigFemScalarResultFrames* calculateStressGradients( int partIndex, const RigFemResultAddress& resVarAddr ); RigFemScalarResultFrames* calculateNodalGradients( int partIndex, const RigFemResultAddress& resVarAddr ); + RigFemScalarResultFrames* calculateNormalizedResult( int partIndex, const RigFemResultAddress& resVarAddr ); const RigFormationNames* activeFormationNames() const; @@ -171,6 +176,7 @@ private: double m_cohesion; double m_frictionAngleRad; + double m_normalizationAirGap; RigStatisticsDataCache* statistics( const RigFemResultAddress& resVarAddr ); std::vector getResAddrToComponentsToRead( const RigFemResultAddress& resVarAddr ); diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemResultAddress.h b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemResultAddress.h index 984a5a9afc..f5afefc7b2 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemResultAddress.h +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemResultAddress.h @@ -35,19 +35,22 @@ public: , componentName( "" ) , timeLapseBaseFrameIdx( NO_TIME_LAPSE ) , refKLayerIndex( NO_COMPACTION ) + , normalizedByHydrostaticPressure( false ) { } RigFemResultAddress( RigFemResultPosEnum resPosType, const std::string& aFieldName, const std::string& aComponentName, - int timeLapseBaseFrameIdx = NO_TIME_LAPSE, - int refKLayerIndex = NO_COMPACTION ) + int timeLapseBaseFrameIdx = NO_TIME_LAPSE, + int refKLayerIndex = NO_COMPACTION, + bool normalizedByHydrostaticPressure = false ) : resultPosType( resPosType ) , fieldName( aFieldName ) , componentName( aComponentName ) , timeLapseBaseFrameIdx( timeLapseBaseFrameIdx ) , refKLayerIndex( refKLayerIndex ) + , normalizedByHydrostaticPressure( normalizedByHydrostaticPressure ) { } @@ -56,6 +59,7 @@ public: std::string componentName; int timeLapseBaseFrameIdx; int refKLayerIndex; + bool normalizedByHydrostaticPressure; static constexpr int allTimeLapsesValue() { @@ -78,6 +82,10 @@ public: { return timeLapseBaseFrameIdx == ALL_TIME_LAPSES; } + bool normalizeByHydrostaticPressure() const + { + return normalizedByHydrostaticPressure; + } bool isValid() const { @@ -92,6 +100,11 @@ public: bool operator<( const RigFemResultAddress& other ) const { + if ( normalizedByHydrostaticPressure != other.normalizedByHydrostaticPressure ) + { + return ( normalizedByHydrostaticPressure < other.normalizedByHydrostaticPressure ); + } + if ( timeLapseBaseFrameIdx != other.timeLapseBaseFrameIdx ) { return ( timeLapseBaseFrameIdx < other.timeLapseBaseFrameIdx ); @@ -118,7 +131,8 @@ public: bool operator==( const RigFemResultAddress& other ) const { if ( resultPosType != other.resultPosType || fieldName != other.fieldName || - componentName != other.componentName || timeLapseBaseFrameIdx != other.timeLapseBaseFrameIdx ) + componentName != other.componentName || timeLapseBaseFrameIdx != other.timeLapseBaseFrameIdx || + normalizedByHydrostaticPressure != other.normalizedByHydrostaticPressure ) { return false; } diff --git a/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.cpp b/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.cpp index 49b874de94..0d49c5c324 100644 --- a/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.cpp +++ b/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.cpp @@ -114,8 +114,8 @@ RimGeoMechResultDefinition::RimGeoMechResultDefinition( void ) "", "", "" ); - CAF_PDM_InitField( &m_normalizationRkbDiff, "NormalizationRkbDiff", -1.0, "Air Gap", "", "", "" ); - m_normalizationRkbDiff.uiCapability()->setUiEditorTypeName( caf::PdmUiDoubleValueEditor::uiEditorTypeName() ); + CAF_PDM_InitField( &m_normalizationAirGap, "NormalizationAirGap", -1.0, "Air Gap", "", "", "" ); + m_normalizationAirGap.uiCapability()->setUiEditorTypeName( caf::PdmUiDoubleValueEditor::uiEditorTypeName() ); CAF_PDM_InitField( &m_compactionRefLayerUiField, "CompactionRefLayerUi", @@ -163,7 +163,7 @@ void RimGeoMechResultDefinition::defineUiOrdering( QString uiConfigName, caf::Pd normalizationGroup->add( &m_normalizeByHydrostaticPressure ); if ( m_normalizeByHydrostaticPressure ) { - normalizationGroup->add( &m_normalizationRkbDiff ); + normalizationGroup->add( &m_normalizationAirGap ); } } @@ -338,30 +338,35 @@ void RimGeoMechResultDefinition::fieldChangedByUi( const caf::PdmFieldHandle* ch m_resultVariableUiField = ""; } } - if ( &m_normalizeByHydrostaticPressure == changedField && m_normalizationRkbDiff < 0.0 ) + if ( &m_normalizeByHydrostaticPressure == changedField && m_normalizationAirGap < 0.0 ) { - RiaMedianCalculator rkbDiffCalc; + RiaMedianCalculator airGapCalc; for ( auto wellPath : RiaApplication::instance()->project()->allWellPaths() ) { if ( wellPath->wellPathGeometry() ) { - double rkbDiff = wellPath->wellPathGeometry()->rkbDiff(); - if ( rkbDiff > 0.0 ) + double airGap = wellPath->wellPathGeometry()->rkbDiff(); + if ( airGap > 0.0 ) { - rkbDiffCalc.add( rkbDiff ); + airGapCalc.add( airGap ); } } } - if ( rkbDiffCalc.valid() ) + if ( airGapCalc.valid() ) { - m_normalizationRkbDiff = rkbDiffCalc.median(); + m_normalizationAirGap = airGapCalc.median(); } else { - m_normalizationRkbDiff = 0.0; + m_normalizationAirGap = 0.0; } } + if ( m_normalizeByHydrostaticPressure && m_normalizationAirGap >= 0.0 ) + { + m_geomCase->geoMechData()->femPartResults()->setNormalizationAirGap( m_normalizationAirGap ); + } + // Get the possible property filter owner RimGeoMechPropertyFilter* propFilter = dynamic_cast( this->parentField()->ownerObject() ); RimGridView* view = nullptr; @@ -372,7 +377,8 @@ void RimGeoMechResultDefinition::fieldChangedByUi( const caf::PdmFieldHandle* ch this->firstAncestorOrThisOfType( rim3dWellLogCurve ); if ( &m_resultVariableUiField == changedField || &m_compactionRefLayerUiField == changedField || - &m_timeLapseBaseTimestep == changedField ) + &m_timeLapseBaseTimestep == changedField || &m_normalizeByHydrostaticPressure == changedField || + &m_normalizationAirGap == changedField ) { QStringList fieldComponentNames = m_resultVariableUiField().split( QRegExp( "\\s+" ) ); if ( fieldComponentNames.size() > 0 ) @@ -535,7 +541,7 @@ void RimGeoMechResultDefinition::defineEditorAttribute( const caf::PdmFieldHandl QString uiConfigName, caf::PdmUiEditorAttribute* attribute ) { - if ( field == &m_normalizationRkbDiff ) + if ( field == &m_normalizationAirGap ) { auto attr = dynamic_cast( attribute ); if ( attr ) @@ -604,7 +610,8 @@ RigFemResultAddress RimGeoMechResultDefinition::resultAddress() const m_timeLapseBaseTimestep(), resultFieldName().toStdString() == RigFemPartResultsCollection::FIELD_NAME_COMPACTION ? m_compactionRefLayer() - : RigFemResultAddress::noCompactionValue() ); + : RigFemResultAddress::noCompactionValue(), + m_normalizeByHydrostaticPressure ); } } @@ -808,18 +815,7 @@ QString RimGeoMechResultDefinition::convertToUiResultFieldName( QString resultFi //-------------------------------------------------------------------------------------------------- bool RimGeoMechResultDefinition::normalizableResultSelected() const { - if ( m_resultPositionType == RIG_NODAL ) - { - return m_resultFieldName == "POR-Bar"; - } - else if ( m_resultPositionType == RIG_ELEMENT_NODAL && ( m_resultFieldName == "ST" || m_resultFieldName == "SE" ) ) - { - static const std::set validComponents = - {"S11", "S22", "S33", "S12", "S13", "S23", "S1", "S2", "S3", "Q", "STM", "SEM"}; - return validComponents.count( m_resultComponentName ); - } - - return false; + return RigFemPartResultsCollection::isNormalizableResult( this->resultAddress() ); } //-------------------------------------------------------------------------------------------------- @@ -827,11 +823,12 @@ bool RimGeoMechResultDefinition::normalizableResultSelected() const //-------------------------------------------------------------------------------------------------- void RimGeoMechResultDefinition::setResultAddress( const RigFemResultAddress& resultAddress ) { - m_resultPositionType = resultAddress.resultPosType; - m_resultFieldName = QString::fromStdString( resultAddress.fieldName ); - m_resultComponentName = QString::fromStdString( resultAddress.componentName ); - m_timeLapseBaseTimestep = resultAddress.timeLapseBaseFrameIdx; - m_compactionRefLayer = resultAddress.refKLayerIndex; + m_resultPositionType = resultAddress.resultPosType; + m_resultFieldName = QString::fromStdString( resultAddress.fieldName ); + m_resultComponentName = QString::fromStdString( resultAddress.componentName ); + m_timeLapseBaseTimestep = resultAddress.timeLapseBaseFrameIdx; + m_compactionRefLayer = resultAddress.refKLayerIndex; + m_normalizeByHydrostaticPressure = resultAddress.normalizedByHydrostaticPressure; m_resultPositionTypeUiField = m_resultPositionType; m_resultVariableUiField = composeFieldCompString( m_resultFieldName(), m_resultComponentName() ); diff --git a/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.h b/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.h index 785ba36b4b..f3489f993d 100644 --- a/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.h +++ b/ApplicationCode/ProjectDataModel/RimGeoMechResultDefinition.h @@ -118,7 +118,7 @@ private: caf::PdmField m_timeLapseBaseTimestep; caf::PdmField m_compactionRefLayer; caf::PdmField m_normalizeByHydrostaticPressure; - caf::PdmField m_normalizationRkbDiff; + caf::PdmField m_normalizationAirGap; // UI Fields only friend class RimGeoMechPropertyFilter; // Property filter needs the ui fields