From 967d8b9bd4962ec97fc2a89016af5d7f1dfaf804 Mon Sep 17 00:00:00 2001 From: Kristian Bendiksen Date: Thu, 30 Jan 2020 09:10:56 +0100 Subject: [PATCH] #4669 Calculate nodal gradients for POR. --- .../RigFemPartResultsCollection.cpp | 106 ++++++++++++++++++ .../RigFemPartResultsCollection.h | 1 + .../GeoMechDataModel/RigHexGradientTools.cpp | 25 ++--- .../GeoMechDataModel/RigHexGradientTools.h | 1 - 4 files changed, 116 insertions(+), 17 deletions(-) diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp index d33c6c0837..6765d357fc 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp @@ -352,6 +352,9 @@ std::map> { fieldCompNames = m_readerInterface->scalarNodeFieldAndComponentNames(); fieldCompNames["POR-Bar"]; + fieldCompNames["PORG"].push_back( "X" ); + fieldCompNames["PORG"].push_back( "Y" ); + fieldCompNames["PORG"].push_back( "Z" ); fieldCompNames[FIELD_NAME_COMPACTION]; } else if ( resPos == RIG_ELEMENT_NODAL ) @@ -1078,6 +1081,103 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateStressGradient( return requestedStress; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemScalarResultFrames* RigFemPartResultsCollection::calculateNodalGradient( int partIndex, + const RigFemResultAddress& resVarAddr ) +{ + CVF_ASSERT( resVarAddr.fieldName == "PORG" ); + CVF_ASSERT( resVarAddr.componentName == "X" || resVarAddr.componentName == "Y" || resVarAddr.componentName == "Z" ); + + caf::ProgressInfo frameCountProgress( this->frameCount() * 5, "" ); + frameCountProgress.setProgressDescription( + "Calculating gradient: " + QString::fromStdString( resVarAddr.fieldName + ": " + resVarAddr.componentName ) ); + frameCountProgress.setNextProgressIncrement( this->frameCount() ); + + RigFemScalarResultFrames* dataFramesX = m_femPartResults[partIndex]->createScalarResult( + RigFemResultAddress( RIG_NODAL, resVarAddr.fieldName, "X" ) ); + frameCountProgress.incrementProgress(); + frameCountProgress.setNextProgressIncrement( this->frameCount() ); + + RigFemScalarResultFrames* dataFramesY = m_femPartResults[partIndex]->createScalarResult( + RigFemResultAddress( RIG_NODAL, resVarAddr.fieldName, "Y" ) ); + frameCountProgress.incrementProgress(); + frameCountProgress.setNextProgressIncrement( this->frameCount() ); + + RigFemScalarResultFrames* dataFramesZ = m_femPartResults[partIndex]->createScalarResult( + RigFemResultAddress( RIG_NODAL, resVarAddr.fieldName, "Z" ) ); + frameCountProgress.incrementProgress(); + frameCountProgress.setNextProgressIncrement( this->frameCount() ); + + RigFemResultAddress porResultAddr( RIG_NODAL, "POR-Bar", "" ); + RigFemScalarResultFrames* srcDataFrames = this->findOrLoadScalarResult( partIndex, porResultAddr ); + + frameCountProgress.incrementProgress(); + + const RigFemPart* femPart = m_femParts->part( partIndex ); + float inf = std::numeric_limits::infinity(); + + 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& dstFrameDataX = dataFramesX->frameData( fIdx ); + std::vector& dstFrameDataY = dataFramesY->frameData( fIdx ); + std::vector& dstFrameDataZ = dataFramesZ->frameData( fIdx ); + + if ( srcFrameData.empty() ) continue; // Create empty results if we have no POR result. + + size_t valCount = femPart->elementNodeResultCount(); + dstFrameDataX.resize( valCount, inf ); + dstFrameDataY.resize( valCount, inf ); + dstFrameDataZ.resize( valCount, inf ); + + int elementCount = femPart->elementCount(); + + for ( int elmIdx = 0; elmIdx < elementCount; ++elmIdx ) + { + RigElementType elmType = femPart->elementType( elmIdx ); + + if ( !( elmType == HEX8P ) ) continue; + + // Find the corner coordinates and values for the node + std::array hexCorners; + std::array cornerValues; + + int elmNodeCount = RigFemTypes::elmentNodeCount( elmType ); + for ( int elmNodIdx = 0; elmNodIdx < elmNodeCount; ++elmNodIdx ) + { + size_t elmNodResIdx = femPart->elementNodeResultIdx( elmIdx, elmNodIdx ); + int nodeIdx = femPart->resultValueIdxFromResultPosType( RIG_NODAL, elmIdx, elmNodIdx ); + + cornerValues[elmNodIdx] = srcFrameData[nodeIdx]; + hexCorners[elmNodIdx] = cvf::Vec3d( nodeCoords[nodeIdx] ); + } + + std::array gradients = RigHexGradientTools::gradients( hexCorners, cornerValues ); + + for ( int elmNodIdx = 0; elmNodIdx < elmNodeCount; ++elmNodIdx ) + { + size_t elmNodResIdx = femPart->elementNodeResultIdx( elmIdx, elmNodIdx ); + int nodeIdx = femPart->resultValueIdxFromResultPosType( RIG_NODAL, elmIdx, elmNodIdx ); + + dstFrameDataX[nodeIdx] = gradients[elmNodIdx].x(); + dstFrameDataY[nodeIdx] = gradients[elmNodIdx].y(); + dstFrameDataZ[nodeIdx] = gradients[elmNodIdx].z(); + } + } + + frameCountProgress.incrementProgress(); + } + + RigFemScalarResultFrames* requestedGradient = this->findOrLoadScalarResult( partIndex, resVarAddr ); + CVF_ASSERT( requestedGradient ); + return requestedGradient; +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -2287,6 +2387,12 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateDerivedResult( i return calculateEnIpPorBarResult( partIndex, resVarAddr ); } + if ( resVarAddr.fieldName == "PORG" && + ( resVarAddr.componentName == "X" || resVarAddr.componentName == "Y" || resVarAddr.componentName == "Z" ) ) + { + return calculateNodalGradient( partIndex, resVarAddr ); + } + if ( ( resVarAddr.fieldName == "NE" ) && ( resVarAddr.componentName == "E11" || resVarAddr.componentName == "E22" || resVarAddr.componentName == "E33" || resVarAddr.componentName == "E12" || resVarAddr.componentName == "E13" || resVarAddr.componentName == "E23" ) ) diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h index 2d81e1ab53..6687b4d927 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h @@ -157,6 +157,7 @@ private: RigFemScalarResultFrames* calculateGamma( int partIndex, const RigFemResultAddress& resVarAddr ); RigFemScalarResultFrames* calculateFormationIndices( int partIndex, const RigFemResultAddress& resVarAddr ); RigFemScalarResultFrames* calculateStressGradient( int partIndex, const RigFemResultAddress& resVarAddr ); + RigFemScalarResultFrames* calculateNodalGradient( int partIndex, const RigFemResultAddress& resVarAddr ); const RigFormationNames* activeFormationNames() const; diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigHexGradientTools.cpp b/ApplicationCode/GeoMech/GeoMechDataModel/RigHexGradientTools.cpp index 00f76232ba..605cb67106 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigHexGradientTools.cpp +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigHexGradientTools.cpp @@ -29,13 +29,13 @@ std::array RigHexGradientTools::gradients( const std::array gradientsUVW; gradientsUVW[0] = cvf::Vec3d( forwardFD( values, 0, 1 ), forwardFD( values, 0, 3 ), forwardFD( values, 0, 4 ) ); - gradientsUVW[1] = cvf::Vec3d( backwardFD( values, 0, 1 ), forwardFD( values, 1, 2 ), forwardFD( values, 1, 5 ) ); - gradientsUVW[2] = cvf::Vec3d( backwardFD( values, 3, 2 ), backwardFD( values, 1, 2 ), forwardFD( values, 2, 6 ) ); - gradientsUVW[3] = cvf::Vec3d( forwardFD( values, 3, 2 ), backwardFD( values, 0, 3 ), forwardFD( values, 3, 7 ) ); - gradientsUVW[4] = cvf::Vec3d( forwardFD( values, 4, 5 ), forwardFD( values, 4, 7 ), backwardFD( values, 0, 4 ) ); - gradientsUVW[5] = cvf::Vec3d( backwardFD( values, 4, 5 ), forwardFD( values, 5, 6 ), backwardFD( values, 1, 5 ) ); - gradientsUVW[6] = cvf::Vec3d( backwardFD( values, 7, 6 ), backwardFD( values, 5, 6 ), backwardFD( values, 2, 6 ) ); - gradientsUVW[7] = cvf::Vec3d( forwardFD( values, 7, 6 ), backwardFD( values, 4, 7 ), backwardFD( values, 3, 7 ) ); + gradientsUVW[1] = cvf::Vec3d( forwardFD( values, 0, 1 ), forwardFD( values, 1, 2 ), forwardFD( values, 1, 5 ) ); + gradientsUVW[2] = cvf::Vec3d( forwardFD( values, 3, 2 ), forwardFD( values, 1, 2 ), forwardFD( values, 2, 6 ) ); + gradientsUVW[3] = cvf::Vec3d( forwardFD( values, 2, 3 ), forwardFD( values, 0, 3 ), forwardFD( values, 3, 7 ) ); + gradientsUVW[4] = cvf::Vec3d( forwardFD( values, 4, 5 ), forwardFD( values, 4, 7 ), forwardFD( values, 0, 4 ) ); + gradientsUVW[5] = cvf::Vec3d( forwardFD( values, 4, 5 ), forwardFD( values, 5, 6 ), forwardFD( values, 1, 5 ) ); + gradientsUVW[6] = cvf::Vec3d( forwardFD( values, 7, 6 ), forwardFD( values, 5, 6 ), forwardFD( values, 2, 6 ) ); + gradientsUVW[7] = cvf::Vec3d( forwardFD( values, 7, 6 ), forwardFD( values, 4, 7 ), forwardFD( values, 3, 7 ) ); std::array NC; NC[0] = {-1, -1, -1}; @@ -53,7 +53,8 @@ std::array RigHexGradientTools::gradients( const std::array& values, int double h = 2.0; return ( values[to] - values[from] ) / h; } - -//-------------------------------------------------------------------------------------------------- -/// -//-------------------------------------------------------------------------------------------------- -double RigHexGradientTools::backwardFD( const std::array& values, int from, int to ) -{ - return forwardFD( values, to, from ); -} diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigHexGradientTools.h b/ApplicationCode/GeoMech/GeoMechDataModel/RigHexGradientTools.h index ed8a38b0b0..b6484c1b71 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigHexGradientTools.h +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigHexGradientTools.h @@ -39,5 +39,4 @@ private: RigHexGradientTools(); static double forwardFD( const std::array& values, int from, int to ); - static double backwardFD( const std::array& values, int from, int to ); };