mirror of
https://github.com/OPM/ResInsight.git
synced 2025-02-25 18:55:39 -06:00
Merge pull request #5454 from OPM/calculate-stress-gradients-4669
Calculate stress gradients 4669
This commit is contained in:
Kristian Bendiksen
GitHub
commit
d52230fda1
@ -33,9 +33,11 @@ add_library( ${PROJECT_NAME}
|
||||
RigFemResultAddress.h
|
||||
RigFemResultPosEnum.h
|
||||
RimFemResultObserver.h
|
||||
RimFemResultObserver.cpp
|
||||
RimGeoMechGeometrySelectionItem.h
|
||||
RimGeoMechGeometrySelectionItem.cpp
|
||||
RimFemResultObserver.cpp
|
||||
RigHexGradientTools.h
|
||||
RigHexGradientTools.cpp
|
||||
RimGeoMechGeometrySelectionItem.h
|
||||
RimGeoMechGeometrySelectionItem.cpp
|
||||
)
|
||||
|
||||
target_include_directories(${PROJECT_NAME}
|
||||
|
||||
@ -36,6 +36,7 @@
|
||||
#include "RigFemPartResults.h"
|
||||
#include "RigFemScalarResultFrames.h"
|
||||
#include "RigFormationNames.h"
|
||||
#include "RigHexGradientTools.h"
|
||||
#include "RigHexIntersectionTools.h"
|
||||
#include "RigStatisticsDataCache.h"
|
||||
#include "RigWbsParameter.h"
|
||||
@ -45,11 +46,10 @@
|
||||
#include "RimWellLogPlot.h"
|
||||
#include "RimWellLogPlotCollection.h"
|
||||
|
||||
#include "cafProgressInfo.h"
|
||||
#include "cvfBoundingBox.h"
|
||||
|
||||
#include "cafProgressInfo.h"
|
||||
#include "cafTensor3.h"
|
||||
|
||||
#include "cvfBoundingBox.h"
|
||||
#include "cvfGeometryTools.h"
|
||||
#include "cvfMath.h"
|
||||
|
||||
@ -346,12 +346,18 @@ std::map<std::string, std::vector<std::string>>
|
||||
if ( activeFormationNames() ) fieldCompNames["Active Formation Names"];
|
||||
}
|
||||
|
||||
const std::vector<std::string> stressComponentNames = getStressComponentNames();
|
||||
const std::vector<std::string> stressGradientComponentNames = getStressGradientComponentNames();
|
||||
|
||||
if ( m_readerInterface.notNull() )
|
||||
{
|
||||
if ( resPos == RIG_NODAL )
|
||||
{
|
||||
fieldCompNames = m_readerInterface->scalarNodeFieldAndComponentNames();
|
||||
fieldCompNames["POR-Bar"];
|
||||
fieldCompNames["POR_GRADIENTS"].push_back( "X" );
|
||||
fieldCompNames["POR_GRADIENTS"].push_back( "Y" );
|
||||
fieldCompNames["POR_GRADIENTS"].push_back( "Z" );
|
||||
fieldCompNames[FIELD_NAME_COMPACTION];
|
||||
}
|
||||
else if ( resPos == RIG_ELEMENT_NODAL )
|
||||
@ -363,15 +369,10 @@ std::map<std::string, std::vector<std::string>>
|
||||
fieldCompNames["SE"].push_back( "DSM" );
|
||||
fieldCompNames["SE"].push_back( "FOS" );
|
||||
|
||||
fieldCompNames["SE"].push_back( "S11" );
|
||||
fieldCompNames["SE"].push_back( "S22" );
|
||||
fieldCompNames["SE"].push_back( "S33" );
|
||||
fieldCompNames["SE"].push_back( "S12" );
|
||||
fieldCompNames["SE"].push_back( "S13" );
|
||||
fieldCompNames["SE"].push_back( "S23" );
|
||||
fieldCompNames["SE"].push_back( "S1" );
|
||||
fieldCompNames["SE"].push_back( "S2" );
|
||||
fieldCompNames["SE"].push_back( "S3" );
|
||||
for ( auto& s : stressComponentNames )
|
||||
{
|
||||
fieldCompNames["SE"].push_back( s );
|
||||
}
|
||||
|
||||
fieldCompNames["SE"].push_back( "S1inc" );
|
||||
fieldCompNames["SE"].push_back( "S1azi" );
|
||||
@ -380,18 +381,18 @@ std::map<std::string, std::vector<std::string>>
|
||||
fieldCompNames["SE"].push_back( "S3inc" );
|
||||
fieldCompNames["SE"].push_back( "S3azi" );
|
||||
|
||||
for ( auto& s : stressGradientComponentNames )
|
||||
{
|
||||
fieldCompNames["SE_GRADIENTS"].push_back( s );
|
||||
}
|
||||
|
||||
fieldCompNames["ST"].push_back( "STM" );
|
||||
fieldCompNames["ST"].push_back( "Q" );
|
||||
|
||||
fieldCompNames["ST"].push_back( "S11" );
|
||||
fieldCompNames["ST"].push_back( "S22" );
|
||||
fieldCompNames["ST"].push_back( "S33" );
|
||||
fieldCompNames["ST"].push_back( "S12" );
|
||||
fieldCompNames["ST"].push_back( "S13" );
|
||||
fieldCompNames["ST"].push_back( "S23" );
|
||||
fieldCompNames["ST"].push_back( "S1" );
|
||||
fieldCompNames["ST"].push_back( "S2" );
|
||||
fieldCompNames["ST"].push_back( "S3" );
|
||||
for ( auto& s : stressComponentNames )
|
||||
{
|
||||
fieldCompNames["ST"].push_back( s );
|
||||
}
|
||||
|
||||
fieldCompNames["ST"].push_back( "S1inc" );
|
||||
fieldCompNames["ST"].push_back( "S1azi" );
|
||||
@ -400,6 +401,11 @@ std::map<std::string, std::vector<std::string>>
|
||||
fieldCompNames["ST"].push_back( "S3inc" );
|
||||
fieldCompNames["ST"].push_back( "S3azi" );
|
||||
|
||||
for ( auto& s : stressGradientComponentNames )
|
||||
{
|
||||
fieldCompNames["ST_GRADIENTS"].push_back( s );
|
||||
}
|
||||
|
||||
fieldCompNames["Gamma"].push_back( "Gamma1" );
|
||||
fieldCompNames["Gamma"].push_back( "Gamma2" );
|
||||
fieldCompNames["Gamma"].push_back( "Gamma3" );
|
||||
@ -979,6 +985,221 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateMeanStressSTM( i
|
||||
return dstDataFrames;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
RigFemScalarResultFrames* RigFemPartResultsCollection::calculateStressGradients( int partIndex,
|
||||
const RigFemResultAddress& resVarAddr )
|
||||
{
|
||||
CVF_ASSERT( resVarAddr.fieldName == "ST_GRADIENTS" || resVarAddr.fieldName == "SE_GRADIENTS" );
|
||||
|
||||
caf::ProgressInfo frameCountProgress( this->frameCount() * 2, "" );
|
||||
frameCountProgress.setProgressDescription(
|
||||
"Calculating gradient: " + QString::fromStdString( resVarAddr.fieldName + ": " + resVarAddr.componentName ) );
|
||||
frameCountProgress.setNextProgressIncrement( this->frameCount() );
|
||||
|
||||
QString origFieldName = QString::fromStdString( resVarAddr.fieldName );
|
||||
// Strip away "_GRADIENTS" to underlying field data
|
||||
QString underlyingFieldName = origFieldName.left( origFieldName.lastIndexOf( QChar( '_' ) ) );
|
||||
QString origComponentName = QString::fromStdString( resVarAddr.componentName );
|
||||
QString componentName = origComponentName.left( origComponentName.lastIndexOf( QChar( '-' ) ) );
|
||||
|
||||
RigFemScalarResultFrames* inputResultFrames =
|
||||
this->findOrLoadScalarResult( partIndex,
|
||||
RigFemResultAddress( resVarAddr.resultPosType,
|
||||
underlyingFieldName.toStdString(),
|
||||
componentName.toStdString() ) );
|
||||
|
||||
RigFemScalarResultFrames* dataFramesX = m_femPartResults[partIndex]->createScalarResult(
|
||||
RigFemResultAddress( resVarAddr.resultPosType, resVarAddr.fieldName, componentName.toStdString() + "-X" ) );
|
||||
RigFemScalarResultFrames* dataFramesY = m_femPartResults[partIndex]->createScalarResult(
|
||||
RigFemResultAddress( resVarAddr.resultPosType, resVarAddr.fieldName, componentName.toStdString() + "-Y" ) );
|
||||
RigFemScalarResultFrames* dataFramesZ = m_femPartResults[partIndex]->createScalarResult(
|
||||
RigFemResultAddress( resVarAddr.resultPosType, resVarAddr.fieldName, componentName.toStdString() + "-Z" ) );
|
||||
frameCountProgress.incrementProgress();
|
||||
|
||||
const RigFemPart* femPart = m_femParts->part( partIndex );
|
||||
int elementCount = femPart->elementCount();
|
||||
const std::vector<cvf::Vec3f>& nodeCoords = femPart->nodes().coordinates;
|
||||
|
||||
int frameCount = inputResultFrames->frameCount();
|
||||
for ( int fIdx = 0; fIdx < frameCount; ++fIdx )
|
||||
{
|
||||
const std::vector<float>& inputData = inputResultFrames->frameData( fIdx );
|
||||
|
||||
std::vector<float>& dstFrameDataX = dataFramesX->frameData( fIdx );
|
||||
std::vector<float>& dstFrameDataY = dataFramesY->frameData( fIdx );
|
||||
std::vector<float>& dstFrameDataZ = dataFramesZ->frameData( fIdx );
|
||||
size_t valCount = inputData.size();
|
||||
dstFrameDataX.resize( valCount );
|
||||
dstFrameDataY.resize( valCount );
|
||||
dstFrameDataZ.resize( valCount );
|
||||
|
||||
#pragma omp parallel for schedule( dynamic )
|
||||
for ( int elmIdx = 0; elmIdx < elementCount; ++elmIdx )
|
||||
{
|
||||
const int* cornerIndices = femPart->connectivities( elmIdx );
|
||||
RigElementType elmType = femPart->elementType( elmIdx );
|
||||
|
||||
if ( !( elmType == HEX8P || elmType == HEX8 ) ) continue;
|
||||
|
||||
// Find the corner coordinates for element
|
||||
std::array<cvf::Vec3d, 8> hexCorners;
|
||||
for ( int i = 0; i < 8; i++ )
|
||||
{
|
||||
hexCorners[i] = cvf::Vec3d( nodeCoords[cornerIndices[i]] );
|
||||
}
|
||||
|
||||
// Find the corresponding corner values for the element
|
||||
std::array<double, 8> cornerValues;
|
||||
|
||||
int elmNodeCount = RigFemTypes::elmentNodeCount( elmType );
|
||||
for ( int elmNodIdx = 0; elmNodIdx < elmNodeCount; ++elmNodIdx )
|
||||
{
|
||||
size_t elmNodResIdx = femPart->elementNodeResultIdx( elmIdx, elmNodIdx );
|
||||
cornerValues[elmNodIdx] = inputData[elmNodResIdx];
|
||||
}
|
||||
|
||||
std::array<cvf::Vec3d, 8> gradients = RigHexGradientTools::gradients( hexCorners, cornerValues );
|
||||
|
||||
for ( int elmNodIdx = 0; elmNodIdx < elmNodeCount; ++elmNodIdx )
|
||||
{
|
||||
size_t elmNodResIdx = femPart->elementNodeResultIdx( elmIdx, elmNodIdx );
|
||||
dstFrameDataX[elmNodResIdx] = gradients[elmNodIdx].x();
|
||||
dstFrameDataY[elmNodResIdx] = gradients[elmNodIdx].y();
|
||||
dstFrameDataZ[elmNodResIdx] = gradients[elmNodIdx].z();
|
||||
}
|
||||
}
|
||||
|
||||
frameCountProgress.incrementProgress();
|
||||
}
|
||||
|
||||
RigFemScalarResultFrames* requestedStress = this->findOrLoadScalarResult( partIndex, resVarAddr );
|
||||
CVF_ASSERT( requestedStress );
|
||||
return requestedStress;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
RigFemScalarResultFrames* RigFemPartResultsCollection::calculateNodalGradients( int partIndex,
|
||||
const RigFemResultAddress& resVarAddr )
|
||||
{
|
||||
CVF_ASSERT( resVarAddr.fieldName == "POR_GRADIENTS" );
|
||||
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;
|
||||
size_t nodeCount = femPart->nodes().nodeIds.size();
|
||||
|
||||
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();
|
||||
|
||||
#pragma omp parallel for schedule( dynamic )
|
||||
for ( long nodeIdx = 0; nodeIdx < static_cast<long>( nodeCount ); nodeIdx++ )
|
||||
{
|
||||
const std::vector<int> elements = femPart->elementsUsingNode( nodeIdx );
|
||||
|
||||
// Compute the average of the elements contributing to this node
|
||||
cvf::Vec3d result = cvf::Vec3d::ZERO;
|
||||
int nValidElements = 0;
|
||||
for ( int elmIdx : elements )
|
||||
{
|
||||
RigElementType elmType = femPart->elementType( elmIdx );
|
||||
if ( elmType == HEX8P )
|
||||
{
|
||||
// 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 );
|
||||
size_t resultValueIdx = femPart->resultValueIdxFromResultPosType( RIG_NODAL, elmIdx, elmNodIdx );
|
||||
|
||||
cornerValues[elmNodIdx] = srcFrameData[resultValueIdx];
|
||||
hexCorners[elmNodIdx] = cvf::Vec3d( nodeCoords[resultValueIdx] );
|
||||
}
|
||||
|
||||
std::array<cvf::Vec3d, 8> gradients = RigHexGradientTools::gradients( hexCorners, cornerValues );
|
||||
|
||||
for ( int elmNodIdx = 0; elmNodIdx < elmNodeCount; ++elmNodIdx )
|
||||
{
|
||||
size_t elmNodResIdx = femPart->elementNodeResultIdx( elmIdx, elmNodIdx );
|
||||
size_t resultValueIdx = femPart->resultValueIdxFromResultPosType( RIG_NODAL, elmIdx, elmNodIdx );
|
||||
// Only use the gradient for particular corner corresponding to the node
|
||||
if ( static_cast<size_t>( nodeIdx ) == resultValueIdx )
|
||||
{
|
||||
result.add( gradients[elmNodIdx] );
|
||||
}
|
||||
}
|
||||
|
||||
nValidElements++;
|
||||
}
|
||||
}
|
||||
|
||||
if ( nValidElements > 0 )
|
||||
{
|
||||
// Round very small values to zero to avoid ugly coloring when gradients
|
||||
// are dominated by floating point math artifacts.
|
||||
double epsilon = 1.0e-6;
|
||||
result /= static_cast<double>( nValidElements );
|
||||
dstFrameDataX[nodeIdx] = std::abs( result.x() ) > epsilon ? result.x() : 0.0;
|
||||
dstFrameDataY[nodeIdx] = std::abs( result.y() ) > epsilon ? result.y() : 0.0;
|
||||
dstFrameDataZ[nodeIdx] = std::abs( result.z() ) > epsilon ? result.z() : 0.0;
|
||||
}
|
||||
}
|
||||
|
||||
frameCountProgress.incrementProgress();
|
||||
}
|
||||
|
||||
RigFemScalarResultFrames* requestedGradient = this->findOrLoadScalarResult( partIndex, resVarAddr );
|
||||
CVF_ASSERT( requestedGradient );
|
||||
return requestedGradient;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
@ -2160,6 +2381,19 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateDerivedResult( i
|
||||
return calculateMeanStressSTM( partIndex, resVarAddr );
|
||||
}
|
||||
|
||||
if ( resVarAddr.fieldName == "ST_GRADIENTS" || resVarAddr.fieldName == "SE_GRADIENTS" )
|
||||
{
|
||||
const std::vector<std::string> allowedComponentNames = getStressGradientComponentNames();
|
||||
|
||||
for ( auto& allowedComponentName : allowedComponentNames )
|
||||
{
|
||||
if ( resVarAddr.componentName == allowedComponentName )
|
||||
{
|
||||
return calculateStressGradients( partIndex, resVarAddr );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if ( resVarAddr.fieldName == "SE" && resVarAddr.componentName == "SEM" )
|
||||
{
|
||||
return calculateMeanStressSEM( partIndex, resVarAddr );
|
||||
@ -2178,6 +2412,12 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateDerivedResult( i
|
||||
return calculateEnIpPorBarResult( partIndex, resVarAddr );
|
||||
}
|
||||
|
||||
if ( resVarAddr.fieldName == "POR_GRADIENTS" &&
|
||||
( resVarAddr.componentName == "X" || resVarAddr.componentName == "Y" || resVarAddr.componentName == "Z" ) )
|
||||
{
|
||||
return calculateNodalGradients( partIndex, resVarAddr );
|
||||
}
|
||||
|
||||
if ( ( resVarAddr.fieldName == "NE" ) && ( resVarAddr.componentName == "E11" || resVarAddr.componentName == "E22" ||
|
||||
resVarAddr.componentName == "E33" || resVarAddr.componentName == "E12" ||
|
||||
resVarAddr.componentName == "E13" || resVarAddr.componentName == "E23" ) )
|
||||
@ -3060,3 +3300,31 @@ void findReferenceElementForNode( const RigFemPart& part, size_t nodeIdx, size_t
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::vector<std::string> RigFemPartResultsCollection::getStressComponentNames()
|
||||
{
|
||||
return {"S11", "S22", "S33", "S12", "S13", "S23", "S1", "S2", "S3"};
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::vector<std::string> RigFemPartResultsCollection::getStressGradientComponentNames()
|
||||
{
|
||||
std::vector<std::string> directions = {"X", "Y", "Z"};
|
||||
std::vector<std::string> stressComponentNames = getStressComponentNames();
|
||||
|
||||
std::vector<std::string> stressGradientComponentNames;
|
||||
for ( auto& s : stressComponentNames )
|
||||
{
|
||||
for ( auto& d : directions )
|
||||
{
|
||||
stressGradientComponentNames.push_back( s + "-" + d );
|
||||
}
|
||||
}
|
||||
|
||||
return stressGradientComponentNames;
|
||||
}
|
||||
|
||||
@ -29,6 +29,7 @@
|
||||
#include <QString>
|
||||
|
||||
#include <map>
|
||||
#include <string>
|
||||
#include <vector>
|
||||
|
||||
class RifGeoMechReaderInterface;
|
||||
@ -156,6 +157,8 @@ private:
|
||||
RigFemScalarResultFrames* calculateST_12_13_23( int partIndex, const RigFemResultAddress& resVarAddr );
|
||||
RigFemScalarResultFrames* calculateGamma( int partIndex, const RigFemResultAddress& resVarAddr );
|
||||
RigFemScalarResultFrames* calculateFormationIndices( int partIndex, const RigFemResultAddress& resVarAddr );
|
||||
RigFemScalarResultFrames* calculateStressGradients( int partIndex, const RigFemResultAddress& resVarAddr );
|
||||
RigFemScalarResultFrames* calculateNodalGradients( int partIndex, const RigFemResultAddress& resVarAddr );
|
||||
|
||||
const RigFormationNames* activeFormationNames() const;
|
||||
|
||||
@ -172,6 +175,9 @@ private:
|
||||
RigStatisticsDataCache* statistics( const RigFemResultAddress& resVarAddr );
|
||||
std::vector<RigFemResultAddress> getResAddrToComponentsToRead( const RigFemResultAddress& resVarAddr );
|
||||
std::map<RigFemResultAddress, cvf::ref<RigStatisticsDataCache>> m_resultStatistics;
|
||||
|
||||
static std::vector<std::string> getStressComponentNames();
|
||||
static std::vector<std::string> getStressGradientComponentNames();
|
||||
};
|
||||
|
||||
class RigFemPart;
|
||||
|
||||
@ -0,0 +1,70 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2020- 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 "RigHexGradientTools.h"
|
||||
|
||||
#include "../cafHexInterpolator/cafHexInterpolator.h"
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
std::array<cvf::Vec3d, 8> RigHexGradientTools::gradients( const std::array<cvf::Vec3d, 8>& hexCorners,
|
||||
const std::array<double, 8>& values )
|
||||
{
|
||||
//
|
||||
std::array<cvf::Vec3d, 8> gradientsUVW;
|
||||
gradientsUVW[0] = cvf::Vec3d( forwardFD( values, 0, 1 ), forwardFD( values, 0, 3 ), forwardFD( values, 0, 4 ) );
|
||||
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<cvf::Vec3d, 8> NC;
|
||||
NC[0] = {-1, -1, -1};
|
||||
NC[1] = {1, -1, -1};
|
||||
NC[2] = {1, 1, -1};
|
||||
NC[3] = {-1, 1, -1};
|
||||
NC[4] = {-1, -1, 1};
|
||||
NC[5] = {1, -1, 1};
|
||||
NC[6] = {1, 1, 1};
|
||||
NC[7] = {-1, 1, 1};
|
||||
|
||||
std::array<cvf::Vec3d, 8> gradientsXYZ;
|
||||
|
||||
for ( int i = 0; i < 8; i++ )
|
||||
{
|
||||
bool isInvertPossible = false;
|
||||
cvf::Mat3d jacobian = caf::HexInterpolator::jacobi( hexCorners, NC[i] );
|
||||
jacobian.transpose();
|
||||
gradientsXYZ[i] = gradientsUVW[i].getTransformedVector( jacobian.getInverted( &isInvertPossible ) );
|
||||
}
|
||||
|
||||
return gradientsXYZ;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
double RigHexGradientTools::forwardFD( const std::array<double, 8>& values, int from, int to )
|
||||
{
|
||||
double h = 2.0;
|
||||
return ( values[to] - values[from] ) / h;
|
||||
}
|
||||
@ -0,0 +1,42 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2020- 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 "RigHexGradientTools.h"
|
||||
|
||||
#include "cvfVector3.h"
|
||||
|
||||
#include <array>
|
||||
|
||||
//==================================================================================================
|
||||
///
|
||||
//==================================================================================================
|
||||
|
||||
class RigHexGradientTools
|
||||
{
|
||||
public:
|
||||
static std::array<cvf::Vec3d, 8> gradients( const std::array<cvf::Vec3d, 8>& hexCorners,
|
||||
const std::array<double, 8>& values );
|
||||
|
||||
private:
|
||||
// Private to avoid instantiation
|
||||
RigHexGradientTools();
|
||||
|
||||
static double forwardFD( const std::array<double, 8>& values, int from, int to );
|
||||
};
|
||||
@ -63,6 +63,7 @@ ${CMAKE_CURRENT_LIST_DIR}/RiaSummaryCurveAnalyzer-Test.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RiaStdStringTools-Test.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RifWellMeasurementReader-Test.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RiaDateStringParser-Test.cpp
|
||||
${CMAKE_CURRENT_LIST_DIR}/RigHexGradientTools-Test.cpp
|
||||
)
|
||||
|
||||
if (RESINSIGHT_ENABLE_GRPC)
|
||||
|
||||
548
ApplicationCode/UnitTests/RigHexGradientTools-Test.cpp
Normal file
548
ApplicationCode/UnitTests/RigHexGradientTools-Test.cpp
Normal file
@ -0,0 +1,548 @@
|
||||
/////////////////////////////////////////////////////////////////////////////////
|
||||
//
|
||||
// Copyright (C) 2011-2012 Statoil ASA, Ceetron AS
|
||||
//
|
||||
// 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 "gtest/gtest.h"
|
||||
|
||||
#include "RigHexGradientTools.h"
|
||||
|
||||
#include "../cafHexInterpolator/cafHexInterpolator.h"
|
||||
|
||||
#include "cvfMatrix4.h"
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
/// GTEST_FILTER="Test_RigHexGradientTools*"
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST( RigHexGradientTools, GradientsForIdentityElement )
|
||||
{
|
||||
// Identity element
|
||||
std::array<cvf::Vec3d, 8> hexCorners = {cvf::Vec3d( -1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, 1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, 1.0 )};
|
||||
|
||||
std::array<double, 8> cornerValues = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
|
||||
|
||||
std::array<cvf::Vec3d, 8> gradients = RigHexGradientTools::gradients( hexCorners, cornerValues );
|
||||
|
||||
for ( int i = 0; i < 8; i++ )
|
||||
{
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[i].x() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[i].y() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[i].z() );
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST( RigHexGradientTools, GradientsForScaledIdentityElement )
|
||||
{
|
||||
// Identity element
|
||||
std::array<cvf::Vec3d, 8> hexCorners = {cvf::Vec3d( -1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, 1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, 1.0 )};
|
||||
|
||||
for ( auto& v : hexCorners )
|
||||
v *= 2.5;
|
||||
|
||||
std::array<double, 8> cornerValues = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
|
||||
|
||||
std::array<cvf::Vec3d, 8> gradients = RigHexGradientTools::gradients( hexCorners, cornerValues );
|
||||
|
||||
for ( int i = 0; i < 8; i++ )
|
||||
{
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[i].x() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[i].y() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[i].z() );
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST( RigHexGradientTools, GradientsForTranslatedIdentityElement )
|
||||
{
|
||||
// Identity element
|
||||
std::array<cvf::Vec3d, 8> hexCorners = {cvf::Vec3d( -1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, 1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, 1.0 )};
|
||||
|
||||
for ( auto& v : hexCorners )
|
||||
v += {3.2, 9.5, -20.3};
|
||||
|
||||
std::array<double, 8> cornerValues = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
|
||||
|
||||
std::array<cvf::Vec3d, 8> gradients = RigHexGradientTools::gradients( hexCorners, cornerValues );
|
||||
|
||||
for ( int i = 0; i < 8; i++ )
|
||||
{
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[i].x() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[i].y() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[i].z() );
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST( RigHexGradientTools, GradientsForRotatedIdentityElement )
|
||||
{
|
||||
// Identity element
|
||||
std::array<cvf::Vec3d, 8> hexCorners = {cvf::Vec3d( -1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, 1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, 1.0 )};
|
||||
|
||||
cvf::Mat4d rot = cvf::Mat4d::fromRotation( {1, 1, 0}, 3.0 );
|
||||
for ( auto& v : hexCorners )
|
||||
v.transformPoint( rot );
|
||||
|
||||
std::array<double, 8> cornerValues = {1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
|
||||
|
||||
std::array<cvf::Vec3d, 8> gradients = RigHexGradientTools::gradients( hexCorners, cornerValues );
|
||||
|
||||
for ( int i = 0; i < 8; i++ )
|
||||
{
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[i].x() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[i].y() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[i].z() );
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST( RigHexGradientTools, GradientsForElement )
|
||||
{
|
||||
// Identity element
|
||||
std::array<cvf::Vec3d, 8> hexCorners = {cvf::Vec3d( -1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, 1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, 1.0 )};
|
||||
|
||||
// Set a higher value in the first corner
|
||||
std::array<double, 8> cornerValues = {2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
|
||||
|
||||
std::array<cvf::Vec3d, 8> gradients = RigHexGradientTools::gradients( hexCorners, cornerValues );
|
||||
|
||||
// The adjacent corners should a non-zero gradient in the direction
|
||||
// towards the first corner
|
||||
EXPECT_DOUBLE_EQ( -0.5, gradients[0].x() );
|
||||
EXPECT_DOUBLE_EQ( -0.5, gradients[0].y() );
|
||||
EXPECT_DOUBLE_EQ( -0.5, gradients[0].z() );
|
||||
|
||||
EXPECT_DOUBLE_EQ( -0.5, gradients[1].x() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[1].y() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[1].z() );
|
||||
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[3].x() );
|
||||
EXPECT_DOUBLE_EQ( -0.5, gradients[3].y() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[3].z() );
|
||||
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[4].x() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[4].y() );
|
||||
EXPECT_DOUBLE_EQ( -0.5, gradients[4].z() );
|
||||
|
||||
// Non-adjacent should be unaffected
|
||||
std::array<int, 4> independentCorners = {2, 5, 7, 6};
|
||||
for ( auto c : independentCorners )
|
||||
{
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[c].x() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[c].y() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[c].z() );
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST( RigHexGradientTools, GradientsForLongElement )
|
||||
{
|
||||
// Identity element
|
||||
std::array<cvf::Vec3d, 8> hexCorners = {cvf::Vec3d( -1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 2.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 2.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 2.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 2.0, 1.0, 1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, 1.0 )};
|
||||
|
||||
// Set a higher value in the first corner
|
||||
std::array<double, 8> cornerValues = {2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
|
||||
|
||||
std::array<cvf::Vec3d, 8> gradients = RigHexGradientTools::gradients( hexCorners, cornerValues );
|
||||
|
||||
// The adjacent corners should a non-zero gradient in the direction
|
||||
// towards the first corner
|
||||
EXPECT_DOUBLE_EQ( -1.0 / 3.0, gradients[0].x() );
|
||||
EXPECT_DOUBLE_EQ( -0.5, gradients[0].y() );
|
||||
EXPECT_DOUBLE_EQ( -0.5, gradients[0].z() );
|
||||
|
||||
EXPECT_DOUBLE_EQ( -1.0 / 3.0, gradients[1].x() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[1].y() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[1].z() );
|
||||
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[3].x() );
|
||||
EXPECT_DOUBLE_EQ( -0.5, gradients[3].y() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[3].z() );
|
||||
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[4].x() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[4].y() );
|
||||
EXPECT_DOUBLE_EQ( -0.5, gradients[4].z() );
|
||||
|
||||
// Non-adjacent should be unaffected
|
||||
std::array<int, 4> independentCorners = {2, 5, 7, 6};
|
||||
for ( auto c : independentCorners )
|
||||
{
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[c].x() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[c].y() );
|
||||
EXPECT_DOUBLE_EQ( 0.0, gradients[c].z() );
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST( RigHexGradientTools, GenerateJacobianForIdentity )
|
||||
{
|
||||
// Identity element
|
||||
std::array<cvf::Vec3d, 8> hexCorners = {cvf::Vec3d( -1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, 1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, 1.0 )};
|
||||
|
||||
for ( int i = 0; i < 8; i++ )
|
||||
{
|
||||
cvf::Mat3d actualJacobian = caf::HexInterpolator::jacobi( hexCorners, hexCorners[i] );
|
||||
|
||||
// Calculate expected jacobian
|
||||
double dx = 2.0;
|
||||
double dy = 2.0;
|
||||
double dz = 2.0;
|
||||
|
||||
double di = 2.0;
|
||||
double dj = 2.0;
|
||||
double dk = 2.0;
|
||||
|
||||
cvf::Mat3d expectedJacobian( dx / di, 0.0, 0.0, 0.0, dy / dj, 0.0, 0.0, 0.0, dz / dk );
|
||||
|
||||
cvf::Mat3d identityMat = cvf::Mat3d::IDENTITY;
|
||||
|
||||
for ( int col = 0; col < 3; col++ )
|
||||
{
|
||||
for ( int row = 0; row < 3; row++ )
|
||||
{
|
||||
EXPECT_DOUBLE_EQ( expectedJacobian.rowCol( row, col ), actualJacobian.rowCol( row, col ) );
|
||||
EXPECT_DOUBLE_EQ( identityMat.rowCol( row, col ), actualJacobian.rowCol( row, col ) );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST( RigHexGradientTools, GenerateJacobianForIdentityScaled )
|
||||
{
|
||||
// Identity element
|
||||
std::array<cvf::Vec3d, 8> normalizedCorners = {cvf::Vec3d( -1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, 1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, 1.0 )};
|
||||
|
||||
double scale = 2.5;
|
||||
|
||||
std::array<cvf::Vec3d, 8> hexCorners;
|
||||
for ( int i = 0; i < 8; i++ )
|
||||
{
|
||||
hexCorners[i] = normalizedCorners[i] * scale;
|
||||
}
|
||||
|
||||
for ( int i = 0; i < 8; i++ )
|
||||
{
|
||||
cvf::Mat3d actualJacobian = caf::HexInterpolator::jacobi( hexCorners, normalizedCorners[i] );
|
||||
|
||||
// Calculate expected jacobian
|
||||
double dx = 5.0;
|
||||
double dy = 5.0;
|
||||
double dz = 5.0;
|
||||
|
||||
double di = 2.0;
|
||||
double dj = 2.0;
|
||||
double dk = 2.0;
|
||||
|
||||
cvf::Mat3d expectedJacobian( dx / di, 0.0, 0.0, 0.0, dy / dj, 0.0, 0.0, 0.0, dz / dk );
|
||||
|
||||
cvf::Mat3d identityMat = cvf::Mat3d::IDENTITY;
|
||||
|
||||
for ( int col = 0; col < 3; col++ )
|
||||
{
|
||||
for ( int row = 0; row < 3; row++ )
|
||||
{
|
||||
EXPECT_DOUBLE_EQ( expectedJacobian.rowCol( row, col ), actualJacobian.rowCol( row, col ) );
|
||||
EXPECT_DOUBLE_EQ( identityMat.rowCol( row, col ) * scale, actualJacobian.rowCol( row, col ) );
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST( RigHexGradientTools, GenerateJacobianForLongElement )
|
||||
{
|
||||
// Try a more complex element
|
||||
std::array<cvf::Vec3d, 8> hexCorners = {cvf::Vec3d( -1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 2.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 2.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 2.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 2.0, 1.0, 1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, 1.0 )};
|
||||
|
||||
cvf::Vec3d corner0( -1.0, -1.0, -1.0 );
|
||||
|
||||
cvf::Mat3d actualJacobian = caf::HexInterpolator::jacobi( hexCorners, corner0 );
|
||||
|
||||
// Calculate expected jacobian
|
||||
double dx = 3.0;
|
||||
double dy = 2.0;
|
||||
double dz = 2.0;
|
||||
|
||||
double di = 2.0;
|
||||
double dj = 2.0;
|
||||
double dk = 2.0;
|
||||
|
||||
cvf::Mat3d expectedJacobian( dx / di, 0.0, 0.0, 0.0, dy / dj, 0.0, 0.0, 0.0, dz / dk );
|
||||
|
||||
for ( int col = 0; col < 3; col++ )
|
||||
{
|
||||
for ( int row = 0; row < 3; row++ )
|
||||
{
|
||||
EXPECT_DOUBLE_EQ( expectedJacobian.rowCol( row, col ), actualJacobian.rowCol( row, col ) );
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST( RigHexGradientTools, GradientsForRotatedElement )
|
||||
{
|
||||
// Try a more complex element
|
||||
std::array<cvf::Vec3d, 8> hexCorners = {cvf::Vec3d( -1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, 1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, 1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, 1.0 )};
|
||||
|
||||
cvf::Vec3d corner0( -1.0, -1.0, -1.0 );
|
||||
|
||||
cvf::Mat4d rot = cvf::Mat4d::fromRotation( {1, 0, 0}, cvf::PI_D / 2.0 );
|
||||
for ( auto& v : hexCorners )
|
||||
v.transformPoint( rot );
|
||||
|
||||
EXPECT_DOUBLE_EQ( 1.0, hexCorners[0].y() );
|
||||
EXPECT_DOUBLE_EQ( -1.0, hexCorners[4].z() );
|
||||
|
||||
// Set a higher value in the first corner
|
||||
std::array<double, 8> cornerValues = {2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
|
||||
|
||||
std::array<cvf::Vec3d, 8> gradients = RigHexGradientTools::gradients( hexCorners, cornerValues );
|
||||
|
||||
EXPECT_NEAR( -0.5, gradients[0].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.5, gradients[0].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( -0.5, gradients[0].z(), 1.0e-10 );
|
||||
|
||||
EXPECT_NEAR( -0.5, gradients[1].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[1].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[1].z(), 1.0e-10 );
|
||||
|
||||
EXPECT_NEAR( 0.0, gradients[2].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[2].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[2].z(), 1.0e-10 );
|
||||
|
||||
EXPECT_NEAR( 0.0, gradients[3].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[3].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( -0.5, gradients[3].z(), 1.0e-10 );
|
||||
|
||||
EXPECT_NEAR( 0.0, gradients[4].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.5, gradients[4].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[4].z(), 1.0e-10 );
|
||||
|
||||
EXPECT_NEAR( 0.0, gradients[5].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[5].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[5].z(), 1.0e-10 );
|
||||
|
||||
EXPECT_NEAR( 0.0, gradients[6].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[6].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[6].z(), 1.0e-10 );
|
||||
|
||||
EXPECT_NEAR( 0.0, gradients[7].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[7].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[7].z(), 1.0e-10 );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST( RigHexGradientTools, GradientsForRotatedAndStretchedElement )
|
||||
{
|
||||
// Try a more complex element
|
||||
std::array<cvf::Vec3d, 8> hexCorners = {cvf::Vec3d( -1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, -1.0, 3.0 ),
|
||||
cvf::Vec3d( 1.0, -1.0, 3.0 ),
|
||||
cvf::Vec3d( 1.0, 1.0, 3.0 ),
|
||||
cvf::Vec3d( -1.0, 1.0, 3.0 )};
|
||||
|
||||
cvf::Vec3d corner0( -1.0, -1.0, -1.0 );
|
||||
|
||||
cvf::Mat4d rot = cvf::Mat4d::fromRotation( {1, 0, 0}, cvf::PI_D / 2.0 );
|
||||
for ( auto& v : hexCorners )
|
||||
v.transformPoint( rot );
|
||||
|
||||
EXPECT_DOUBLE_EQ( 1.0, hexCorners[0].y() );
|
||||
EXPECT_DOUBLE_EQ( -1.0, hexCorners[4].z() );
|
||||
|
||||
// Set a higher value in the first corner
|
||||
std::array<double, 8> cornerValues = {2.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
|
||||
|
||||
std::array<cvf::Vec3d, 8> gradients = RigHexGradientTools::gradients( hexCorners, cornerValues );
|
||||
|
||||
EXPECT_NEAR( -0.5, gradients[0].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.25, gradients[0].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( -0.5, gradients[0].z(), 1.0e-10 );
|
||||
|
||||
EXPECT_NEAR( -0.5, gradients[1].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[1].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[1].z(), 1.0e-10 );
|
||||
|
||||
EXPECT_NEAR( 0.0, gradients[2].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[2].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[2].z(), 1.0e-10 );
|
||||
|
||||
EXPECT_NEAR( 0.0, gradients[3].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[3].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( -0.5, gradients[3].z(), 1.0e-10 );
|
||||
|
||||
EXPECT_NEAR( 0.0, gradients[4].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.25, gradients[4].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[4].z(), 1.0e-10 );
|
||||
|
||||
EXPECT_NEAR( 0.0, gradients[5].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[5].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[5].z(), 1.0e-10 );
|
||||
|
||||
EXPECT_NEAR( 0.0, gradients[6].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[6].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[6].z(), 1.0e-10 );
|
||||
|
||||
EXPECT_NEAR( 0.0, gradients[7].x(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[7].y(), 1.0e-10 );
|
||||
EXPECT_NEAR( 0.0, gradients[7].z(), 1.0e-10 );
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
///
|
||||
//--------------------------------------------------------------------------------------------------
|
||||
TEST( RigHexGradientTools, GenerateJacobianForRotatedElement )
|
||||
{
|
||||
// Try a more complex element
|
||||
std::array<cvf::Vec3d, 8> hexCorners = {cvf::Vec3d( -1.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 2.0, -1.0, -1.0 ),
|
||||
cvf::Vec3d( 2.0, 3.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, 3.0, -1.0 ),
|
||||
cvf::Vec3d( -1.0, -1.0, 4.0 ),
|
||||
cvf::Vec3d( 2.0, -1.0, 4.0 ),
|
||||
cvf::Vec3d( 2.0, 3.0, 4.0 ),
|
||||
cvf::Vec3d( -1.0, 3.0, 4.0 )};
|
||||
|
||||
cvf::Vec3d corner0( -1.0, -1.0, -1.0 );
|
||||
|
||||
cvf::Mat4d rot = cvf::Mat4d::fromRotation( {1, 0, 0}, cvf::PI_D / 2.0 );
|
||||
for ( auto& v : hexCorners )
|
||||
v.transformPoint( rot );
|
||||
|
||||
cvf::Mat3d actualJacobian = caf::HexInterpolator::jacobi( hexCorners, corner0 );
|
||||
|
||||
// Calculate expected jacobian
|
||||
cvf::Vec3d dx_rot = ( hexCorners[1] - hexCorners[0] );
|
||||
cvf::Vec3d dy_rot = ( hexCorners[3] - hexCorners[0] );
|
||||
cvf::Vec3d dz_rot = ( hexCorners[4] - hexCorners[0] );
|
||||
|
||||
double di = 2.0;
|
||||
double dj = 2.0;
|
||||
double dk = 2.0;
|
||||
|
||||
cvf::Mat3d expectedJacobian( dx_rot.x() / di,
|
||||
dx_rot.y() / di,
|
||||
dx_rot.z() / di,
|
||||
dy_rot.x() / dj,
|
||||
dy_rot.y() / dj,
|
||||
dy_rot.z() / dj,
|
||||
dz_rot.x() / dk,
|
||||
dz_rot.y() / dk,
|
||||
dz_rot.z() / dk );
|
||||
expectedJacobian.transpose();
|
||||
|
||||
for ( int row = 0; row < 3; row++ )
|
||||
{
|
||||
for ( int col = 0; col < 3; col++ )
|
||||
{
|
||||
EXPECT_DOUBLE_EQ( expectedJacobian.rowCol( row, col ), actualJacobian.rowCol( row, col ) );
|
||||
}
|
||||
}
|
||||
}
|
||||
@ -220,6 +220,85 @@ public:
|
||||
return interpolationCoeffs(pointInNormElm);
|
||||
}
|
||||
|
||||
|
||||
static cvf::Mat3d jacobi(const std::array<cvf::Vec3d, 8>& hexCorners, const cvf::Vec3d & pointInNormalizedElement)
|
||||
{
|
||||
const double k_1_8 = 0.125;
|
||||
|
||||
double i = pointInNormalizedElement[0];
|
||||
double j = pointInNormalizedElement[1];
|
||||
double k = pointInNormalizedElement[2];
|
||||
|
||||
double ij = i*j;
|
||||
double ik = i*k;
|
||||
double jk = j*k;
|
||||
|
||||
double C0_x = hexCorners[0][0]; double C0_y = hexCorners[0][1]; double C0_z = hexCorners[0][2];
|
||||
double C1_x = hexCorners[1][0]; double C1_y = hexCorners[1][1]; double C1_z = hexCorners[1][2];
|
||||
double C2_x = hexCorners[2][0]; double C2_y = hexCorners[2][1]; double C2_z = hexCorners[2][2];
|
||||
double C3_x = hexCorners[3][0]; double C3_y = hexCorners[3][1]; double C3_z = hexCorners[3][2];
|
||||
double C4_x = hexCorners[4][0]; double C4_y = hexCorners[4][1]; double C4_z = hexCorners[4][2];
|
||||
double C5_x = hexCorners[5][0]; double C5_y = hexCorners[5][1]; double C5_z = hexCorners[5][2];
|
||||
double C6_x = hexCorners[6][0]; double C6_y = hexCorners[6][1]; double C6_z = hexCorners[6][2];
|
||||
double C7_x = hexCorners[7][0]; double C7_y = hexCorners[7][1]; double C7_z = hexCorners[7][2];
|
||||
|
||||
double dN0_di = (- 1 + j + k - jk);
|
||||
double dN1_di = (+ 1 - j - k + jk);
|
||||
double dN2_di = (+ 1 + j - k - jk);
|
||||
double dN3_di = (- 1 - j + k + jk);
|
||||
double dN4_di = (- 1 + j - k + jk);
|
||||
double dN5_di = (+ 1 - j + k - jk);
|
||||
double dN6_di = (+ 1 + j + k + jk);
|
||||
double dN7_di = (- 1 - j - k - jk);
|
||||
|
||||
double dN0_dj = (- 1 + i + k - ik);
|
||||
double dN1_dj = (- 1 - i + k + ik);
|
||||
double dN2_dj = (+ 1 + i - k - ik);
|
||||
double dN3_dj = (+ 1 - i - k + ik);
|
||||
double dN4_dj = (- 1 + i - k + ik);
|
||||
double dN5_dj = (- 1 - i - k - ik);
|
||||
double dN6_dj = (+ 1 + i + k + ik);
|
||||
double dN7_dj = (+ 1 - i + k - ik);
|
||||
|
||||
double dN0_dk = (- 1 + i + j - ij);
|
||||
double dN1_dk = (- 1 - i + j + ij);
|
||||
double dN2_dk = (- 1 - i - j - ij);
|
||||
double dN3_dk = (- 1 + i - j + ij);
|
||||
double dN4_dk = (+ 1 - i - j + ij);
|
||||
double dN5_dk = (+ 1 + i - j - ij);
|
||||
double dN6_dk = (+ 1 + i + j + ij);
|
||||
double dN7_dk = (+ 1 - i + j - ij);
|
||||
|
||||
double dx_di = k_1_8 * ( (dN0_di) * C0_x + (dN1_di) * C1_x + (dN2_di) * C2_x + (dN3_di) * C3_x + (dN4_di) * C4_x + (dN5_di) * C5_x + (dN6_di) * C6_x + (dN7_di) * C7_x );
|
||||
double dx_dj = k_1_8 * ( (dN0_dj) * C0_x + (dN1_dj) * C1_x + (dN2_dj) * C2_x + (dN3_dj) * C3_x + (dN4_dj) * C4_x + (dN5_dj) * C5_x + (dN6_dj) * C6_x + (dN7_dj) * C7_x );
|
||||
double dx_dk = k_1_8 * ( (dN0_dk) * C0_x + (dN1_dk) * C1_x + (dN2_dk) * C2_x + (dN3_dk) * C3_x + (dN4_dk) * C4_x + (dN5_dk) * C5_x + (dN6_dk) * C6_x + (dN7_dk) * C7_x );
|
||||
|
||||
double dy_di = k_1_8 * ( (dN0_di) * C0_y + (dN1_di) * C1_y + (dN2_di) * C2_y + (dN3_di) * C3_y + (dN4_di) * C4_y + (dN5_di) * C5_y + (dN6_di) * C6_y + (dN7_di) * C7_y );
|
||||
double dy_dj = k_1_8 * ( (dN0_dj) * C0_y + (dN1_dj) * C1_y + (dN2_dj) * C2_y + (dN3_dj) * C3_y + (dN4_dj) * C4_y + (dN5_dj) * C5_y + (dN6_dj) * C6_y + (dN7_dj) * C7_y );
|
||||
double dy_dk = k_1_8 * ( (dN0_dk) * C0_y + (dN1_dk) * C1_y + (dN2_dk) * C2_y + (dN3_dk) * C3_y + (dN4_dk) * C4_y + (dN5_dk) * C5_y + (dN6_dk) * C6_y + (dN7_dk) * C7_y );
|
||||
|
||||
double dz_di = k_1_8 * ( (dN0_di) * C0_z + (dN1_di) * C1_z + (dN2_di) * C2_z + (dN3_di) * C3_z + (dN4_di) * C4_z + (dN5_di) * C5_z + (dN6_di) * C6_z + (dN7_di) * C7_z );
|
||||
double dz_dj = k_1_8 * ( (dN0_dj) * C0_z + (dN1_dj) * C1_z + (dN2_dj) * C2_z + (dN3_dj) * C3_z + (dN4_dj) * C4_z + (dN5_dj) * C5_z + (dN6_dj) * C6_z + (dN7_dj) * C7_z );
|
||||
double dz_dk = k_1_8 * ( (dN0_dk) * C0_z + (dN1_dk) * C1_z + (dN2_dk) * C2_z + (dN3_dk) * C3_z + (dN4_dk) * C4_z + (dN5_dk) * C5_z + (dN6_dk) * C6_z + (dN7_dk) * C7_z );
|
||||
|
||||
// Do not know which ordering ends up as correct
|
||||
|
||||
#if 1 // Use literature jacobi ordering
|
||||
return cvf::Mat3d(
|
||||
dx_di, dx_dj, dx_dk,
|
||||
dy_di, dy_dj, dy_dk,
|
||||
dz_di, dz_dj, dz_dk
|
||||
);
|
||||
#else // use NTNU ordering
|
||||
return cvf::Mat3d(
|
||||
dx_di, dy_di, dz_di,
|
||||
dx_dj, dy_dj, dz_dj,
|
||||
dx_dk, dy_dk, dz_dk
|
||||
);
|
||||
#endif
|
||||
}
|
||||
|
||||
|
||||
private:
|
||||
|
||||
static double interpolateInNormElm( const cvf::Vec3d & pointInNormElm, const std::array<double, 8>& values)
|
||||
@ -300,83 +379,6 @@ private:
|
||||
|
||||
return normPoint;
|
||||
}
|
||||
|
||||
static cvf::Mat3d jacobi(const std::array<cvf::Vec3d, 8>& hexCorners, const cvf::Vec3d & pointInNormalizedElement)
|
||||
{
|
||||
const double k_1_8 = 0.125;
|
||||
|
||||
double i = pointInNormalizedElement[0];
|
||||
double j = pointInNormalizedElement[1];
|
||||
double k = pointInNormalizedElement[2];
|
||||
|
||||
double ij = i*j;
|
||||
double ik = i*k;
|
||||
double jk = j*k;
|
||||
|
||||
double C0_x = hexCorners[0][0]; double C0_y = hexCorners[0][1]; double C0_z = hexCorners[0][2];
|
||||
double C1_x = hexCorners[1][0]; double C1_y = hexCorners[1][1]; double C1_z = hexCorners[1][2];
|
||||
double C2_x = hexCorners[2][0]; double C2_y = hexCorners[2][1]; double C2_z = hexCorners[2][2];
|
||||
double C3_x = hexCorners[3][0]; double C3_y = hexCorners[3][1]; double C3_z = hexCorners[3][2];
|
||||
double C4_x = hexCorners[4][0]; double C4_y = hexCorners[4][1]; double C4_z = hexCorners[4][2];
|
||||
double C5_x = hexCorners[5][0]; double C5_y = hexCorners[5][1]; double C5_z = hexCorners[5][2];
|
||||
double C6_x = hexCorners[6][0]; double C6_y = hexCorners[6][1]; double C6_z = hexCorners[6][2];
|
||||
double C7_x = hexCorners[7][0]; double C7_y = hexCorners[7][1]; double C7_z = hexCorners[7][2];
|
||||
|
||||
double dN0_di = (- 1 + j + k - jk);
|
||||
double dN1_di = (+ 1 - j - k + jk);
|
||||
double dN2_di = (+ 1 + j - k - jk);
|
||||
double dN3_di = (- 1 - j + k + jk);
|
||||
double dN4_di = (- 1 + j - k + jk);
|
||||
double dN5_di = (+ 1 - j + k - jk);
|
||||
double dN6_di = (+ 1 + j + k + jk);
|
||||
double dN7_di = (- 1 - j - k - jk);
|
||||
|
||||
double dN0_dj = (- 1 + i + k - ik);
|
||||
double dN1_dj = (- 1 - i + k + ik);
|
||||
double dN2_dj = (+ 1 + i - k - ik);
|
||||
double dN3_dj = (+ 1 - i - k + ik);
|
||||
double dN4_dj = (- 1 + i - k + ik);
|
||||
double dN5_dj = (- 1 - i - k - ik);
|
||||
double dN6_dj = (+ 1 + i + k + ik);
|
||||
double dN7_dj = (+ 1 - i + k - ik);
|
||||
|
||||
double dN0_dk = (- 1 + i + j - ij);
|
||||
double dN1_dk = (- 1 - i + j + ij);
|
||||
double dN2_dk = (- 1 - i - j - ij);
|
||||
double dN3_dk = (- 1 + i - j + ij);
|
||||
double dN4_dk = (+ 1 - i - j + ij);
|
||||
double dN5_dk = (+ 1 + i - j - ij);
|
||||
double dN6_dk = (+ 1 + i + j + ij);
|
||||
double dN7_dk = (+ 1 - i + j - ij);
|
||||
|
||||
double dx_di = k_1_8 * ( (dN0_di) * C0_x + (dN1_di) * C1_x + (dN2_di) * C2_x + (dN3_di) * C3_x + (dN4_di) * C4_x + (dN5_di) * C5_x + (dN6_di) * C6_x + (dN7_di) * C7_x );
|
||||
double dx_dj = k_1_8 * ( (dN0_dj) * C0_x + (dN1_dj) * C1_x + (dN2_dj) * C2_x + (dN3_dj) * C3_x + (dN4_dj) * C4_x + (dN5_dj) * C5_x + (dN6_dj) * C6_x + (dN7_dj) * C7_x );
|
||||
double dx_dk = k_1_8 * ( (dN0_dk) * C0_x + (dN1_dk) * C1_x + (dN2_dk) * C2_x + (dN3_dk) * C3_x + (dN4_dk) * C4_x + (dN5_dk) * C5_x + (dN6_dk) * C6_x + (dN7_dk) * C7_x );
|
||||
|
||||
double dy_di = k_1_8 * ( (dN0_di) * C0_y + (dN1_di) * C1_y + (dN2_di) * C2_y + (dN3_di) * C3_y + (dN4_di) * C4_y + (dN5_di) * C5_y + (dN6_di) * C6_y + (dN7_di) * C7_y );
|
||||
double dy_dj = k_1_8 * ( (dN0_dj) * C0_y + (dN1_dj) * C1_y + (dN2_dj) * C2_y + (dN3_dj) * C3_y + (dN4_dj) * C4_y + (dN5_dj) * C5_y + (dN6_dj) * C6_y + (dN7_dj) * C7_y );
|
||||
double dy_dk = k_1_8 * ( (dN0_dk) * C0_y + (dN1_dk) * C1_y + (dN2_dk) * C2_y + (dN3_dk) * C3_y + (dN4_dk) * C4_y + (dN5_dk) * C5_y + (dN6_dk) * C6_y + (dN7_dk) * C7_y );
|
||||
|
||||
double dz_di = k_1_8 * ( (dN0_di) * C0_z + (dN1_di) * C1_z + (dN2_di) * C2_z + (dN3_di) * C3_z + (dN4_di) * C4_z + (dN5_di) * C5_z + (dN6_di) * C6_z + (dN7_di) * C7_z );
|
||||
double dz_dj = k_1_8 * ( (dN0_dj) * C0_z + (dN1_dj) * C1_z + (dN2_dj) * C2_z + (dN3_dj) * C3_z + (dN4_dj) * C4_z + (dN5_dj) * C5_z + (dN6_dj) * C6_z + (dN7_dj) * C7_z );
|
||||
double dz_dk = k_1_8 * ( (dN0_dk) * C0_z + (dN1_dk) * C1_z + (dN2_dk) * C2_z + (dN3_dk) * C3_z + (dN4_dk) * C4_z + (dN5_dk) * C5_z + (dN6_dk) * C6_z + (dN7_dk) * C7_z );
|
||||
|
||||
// Do not know which ordering ends up as correct
|
||||
|
||||
#if 1 // Use literature jacobi ordering
|
||||
return cvf::Mat3d(
|
||||
dx_di, dx_dj, dx_dk,
|
||||
dy_di, dy_dj, dy_dk,
|
||||
dz_di, dz_dj, dz_dk
|
||||
);
|
||||
#else // use NTNU ordering
|
||||
return cvf::Mat3d(
|
||||
dx_di, dy_di, dz_di,
|
||||
dx_dj, dy_dj, dz_dj,
|
||||
dx_dk, dy_dk, dz_dk
|
||||
);
|
||||
#endif
|
||||
}
|
||||
};
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
Loading…
Reference in New Issue
Block a user