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

#4669 Calculate nodal gradients for POR.

Browse Source
This commit is contained in:
Kristian Bendiksen 2020-01-30 09:10:56 +01:00
parent 5556fb39b1
commit 967d8b9bd4
4 changed files with 116 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

106
ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp
View File

@ -352,6 +352,9 @@ std::map<std::string, std::vector<std::string>>
{ {
fieldCompNames = m_readerInterface->scalarNodeFieldAndComponentNames(); fieldCompNames = m_readerInterface->scalarNodeFieldAndComponentNames();
fieldCompNames["POR-Bar"]; fieldCompNames["POR-Bar"];
fieldCompNames["PORG"].push_back( "X" );
fieldCompNames["PORG"].push_back( "Y" );
fieldCompNames["PORG"].push_back( "Z" );
fieldCompNames[FIELD_NAME_COMPACTION]; fieldCompNames[FIELD_NAME_COMPACTION];
} }
else if ( resPos == RIG_ELEMENT_NODAL ) else if ( resPos == RIG_ELEMENT_NODAL )
@ -1078,6 +1081,103 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateStressGradient(
return requestedStress; 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<float>::infinity();
const std::vector<cvf::Vec3f>& nodeCoords = femPart->nodes().coordinates;
int frameCount = srcDataFrames->frameCount();
for ( int fIdx = 0; fIdx < frameCount; ++fIdx )
{
const std::vector<float>& srcFrameData = srcDataFrames->frameData( fIdx );
std::vector<float>& dstFrameDataX = dataFramesX->frameData( fIdx );
std::vector<float>& dstFrameDataY = dataFramesY->frameData( fIdx );
std::vector<float>& 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<cvf::Vec3d, 8> hexCorners;
std::array<double, 8> 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<cvf::Vec3d, 8> 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 ); 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" || if ( ( resVarAddr.fieldName == "NE" ) && ( resVarAddr.componentName == "E11" || resVarAddr.componentName == "E22" ||
resVarAddr.componentName == "E33" || resVarAddr.componentName == "E12" || resVarAddr.componentName == "E33" || resVarAddr.componentName == "E12" ||
resVarAddr.componentName == "E13" || resVarAddr.componentName == "E23" ) ) resVarAddr.componentName == "E13" || resVarAddr.componentName == "E23" ) )

1
ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h
View File

@ -157,6 +157,7 @@ private:
RigFemScalarResultFrames* calculateGamma( int partIndex, const RigFemResultAddress& resVarAddr ); RigFemScalarResultFrames* calculateGamma( int partIndex, const RigFemResultAddress& resVarAddr );
RigFemScalarResultFrames* calculateFormationIndices( int partIndex, const RigFemResultAddress& resVarAddr ); RigFemScalarResultFrames* calculateFormationIndices( int partIndex, const RigFemResultAddress& resVarAddr );
RigFemScalarResultFrames* calculateStressGradient( int partIndex, const RigFemResultAddress& resVarAddr ); RigFemScalarResultFrames* calculateStressGradient( int partIndex, const RigFemResultAddress& resVarAddr );
RigFemScalarResultFrames* calculateNodalGradient( int partIndex, const RigFemResultAddress& resVarAddr );
const RigFormationNames* activeFormationNames() const; const RigFormationNames* activeFormationNames() const;

25
ApplicationCode/GeoMech/GeoMechDataModel/RigHexGradientTools.cpp
View File

@ -29,13 +29,13 @@ std::array<cvf::Vec3d, 8> RigHexGradientTools::gradients( const std::array<cvf::
// //
std::array<cvf::Vec3d, 8> gradientsUVW; std::array<cvf::Vec3d, 8> gradientsUVW;
gradientsUVW[0] = cvf::Vec3d( forwardFD( values, 0, 1 ), forwardFD( values, 0, 3 ), forwardFD( values, 0, 4 ) ); 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[1] = cvf::Vec3d( forwardFD( 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[2] = cvf::Vec3d( forwardFD( values, 3, 2 ), forwardFD( values, 1, 2 ), forwardFD( values, 2, 6 ) );
gradientsUVW[3] = cvf::Vec3d( forwardFD( values, 3, 2 ), backwardFD( values, 0, 3 ), forwardFD( values, 3, 7 ) ); 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 ), backwardFD( values, 0, 4 ) ); gradientsUVW[4] = cvf::Vec3d( forwardFD( values, 4, 5 ), forwardFD( values, 4, 7 ), forwardFD( values, 0, 4 ) );
gradientsUVW[5] = cvf::Vec3d( backwardFD( values, 4, 5 ), forwardFD( values, 5, 6 ), backwardFD( values, 1, 5 ) ); gradientsUVW[5] = cvf::Vec3d( forwardFD( values, 4, 5 ), forwardFD( values, 5, 6 ), forwardFD( values, 1, 5 ) );
gradientsUVW[6] = cvf::Vec3d( backwardFD( values, 7, 6 ), backwardFD( values, 5, 6 ), backwardFD( values, 2, 6 ) ); gradientsUVW[6] = cvf::Vec3d( forwardFD( values, 7, 6 ), forwardFD( values, 5, 6 ), forwardFD( values, 2, 6 ) );
gradientsUVW[7] = cvf::Vec3d( forwardFD( values, 7, 6 ), backwardFD( values, 4, 7 ), backwardFD( values, 3, 7 ) ); gradientsUVW[7] = cvf::Vec3d( forwardFD( values, 7, 6 ), forwardFD( values, 4, 7 ), forwardFD( values, 3, 7 ) );
std::array<cvf::Vec3d, 8> NC; std::array<cvf::Vec3d, 8> NC;
NC[0] = {-1, -1, -1}; NC[0] = {-1, -1, -1};
@ -53,7 +53,8 @@ std::array<cvf::Vec3d, 8> RigHexGradientTools::gradients( const std::array<cvf::
{ {
bool isInvertPossible = false; bool isInvertPossible = false;
cvf::Mat3d jacobian = caf::HexInterpolator::jacobi( hexCorners, NC[i] ); cvf::Mat3d jacobian = caf::HexInterpolator::jacobi( hexCorners, NC[i] );
gradientsXYZ[i] = gradientsUVW[i].getTransformedVector( jacobian.getInverted( &isInvertPossible ) ); // jacobian.transpose();
gradientsXYZ[i] = gradientsUVW[i].getTransformedVector( jacobian.getInverted( &isInvertPossible ) );
} }
return gradientsXYZ; return gradientsXYZ;
@ -67,11 +68,3 @@ double RigHexGradientTools::forwardFD( const std::array<double, 8>& values, int
double h = 2.0; double h = 2.0;
return ( values[to] - values[from] ) / h; return ( values[to] - values[from] ) / h;
} }
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
double RigHexGradientTools::backwardFD( const std::array<double, 8>& values, int from, int to )
{
return forwardFD( values, to, from );
}

1
ApplicationCode/GeoMech/GeoMechDataModel/RigHexGradientTools.h
View File

@ -39,5 +39,4 @@ private:
RigHexGradientTools(); RigHexGradientTools();
static double forwardFD( const std::array<double, 8>& values, int from, int to ); static double forwardFD( const std::array<double, 8>& values, int from, int to );
static double backwardFD( const std::array<double, 8>& values, int from, int to );
}; };