#941 Added TP, FAULTMOB, PCRIT, TPinc as transformed stress results

This commit is contained in:
Jacob Støren 2016-11-01 09:29:52 +01:00
parent 055158aeea
commit 8d2d4dfe30
4 changed files with 160 additions and 40 deletions

View File

@ -343,16 +343,25 @@ std::map<std::string, std::vector<std::string> > RigFemPartResultsCollection::sc
fieldCompNames["SE"].push_back("SN");
fieldCompNames["SE"].push_back("STH");
fieldCompNames["SE"].push_back("STQV");
fieldCompNames["SE"].push_back("TP");
fieldCompNames["SE"].push_back("TPinc");
fieldCompNames["SE"].push_back("TNH" );
fieldCompNames["SE"].push_back("TNQV");
fieldCompNames["SE"].push_back("THQV");
fieldCompNames["SE"].push_back("FAULTMOB");
fieldCompNames["SE"].push_back("PCRIT");
fieldCompNames["ST"].push_back("SN");
fieldCompNames["ST"].push_back("STH");
fieldCompNames["ST"].push_back("STQV");
fieldCompNames["ST"].push_back("TP");
fieldCompNames["ST"].push_back("TPinc");
fieldCompNames["ST"].push_back("TNH");
fieldCompNames["ST"].push_back("TNQV");
fieldCompNames["ST"].push_back("THQV");
fieldCompNames["ST"].push_back("FAULTMOB");
fieldCompNames["ST"].push_back("PCRIT");
}
}
@ -646,7 +655,7 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateDSM(int partInde
frameCountProgress.incrementProgress();
float tanFricAng = tan(m_frictionAngleRad);
float cohPrFricAngle = (float)(m_cohesion/tanFricAng);
float cohPrTanFricAngle = (float)(m_cohesion/tanFricAng);
int frameCount = se1Frames->frameCount();
for ( int fIdx = 0; fIdx < frameCount; ++fIdx )
{
@ -661,7 +670,7 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateDSM(int partInde
{
float se1 = se1Data[vIdx];
float se3 = se3Data[vIdx];
float rho = 2.0f * atan( sqrt(( se1 + cohPrFricAngle)/(se3 + cohPrFricAngle)) - M_PI_4);
float rho = 2.0f * atan( sqrt(( se1 + cohPrTanFricAngle)/(se3 + cohPrTanFricAngle)) - M_PI_4);
{
dstFrameData[vIdx] = tan(rho)/tanFricAng;
@ -915,7 +924,9 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateDeviatoricStrain
RigFemScalarResultFrames* RigFemPartResultsCollection::calculateSurfaceAlignedStress(int partIndex, const RigFemResultAddress& resVarAddr)
{
CVF_ASSERT( resVarAddr.componentName == "STH" || resVarAddr.componentName == "STQV" || resVarAddr.componentName == "SN"
|| resVarAddr.componentName == "TNH" || resVarAddr.componentName == "TNQV" || resVarAddr.componentName == "THQV");
|| resVarAddr.componentName == "TNH" || resVarAddr.componentName == "TNQV" || resVarAddr.componentName == "THQV"
|| resVarAddr.componentName == "TP" || resVarAddr.componentName == "TPinc"
|| resVarAddr.componentName == "FAULTMOB" || resVarAddr.componentName == "PCRIT");
caf::ProgressInfo frameCountProgress(this->frameCount() * 7, "");
frameCountProgress.setProgressDescription("Calculating " + QString::fromStdString(resVarAddr.fieldName + ": " + resVarAddr.componentName));
@ -933,18 +944,25 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateSurfaceAlignedSt
frameCountProgress.incrementProgress(); frameCountProgress.setNextProgressIncrement(this->frameCount());
RigFemScalarResultFrames * s13Frames = this->findOrLoadScalarResult(partIndex, RigFemResultAddress(RIG_ELEMENT_NODAL, resVarAddr.fieldName, "S13"));
RigFemScalarResultFrames * sNormFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "SN"));
RigFemScalarResultFrames * sHoriFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "STH"));
RigFemScalarResultFrames * sVertFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "STQV"));
RigFemScalarResultFrames * sNHFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "TNH" ));
RigFemScalarResultFrames * sNVFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "TNQV"));
RigFemScalarResultFrames * sHVFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "THQV"));
RigFemScalarResultFrames * SNFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "SN"));
RigFemScalarResultFrames * STHFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "STH"));
RigFemScalarResultFrames * STQVFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "STQV"));
RigFemScalarResultFrames * TNHFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "TNH" ));
RigFemScalarResultFrames * TNQVFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "TNQV"));
RigFemScalarResultFrames * THQVFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "THQV"));
RigFemScalarResultFrames * TPFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "TP"));
RigFemScalarResultFrames * TPincFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "TPinc"));
RigFemScalarResultFrames * FAULTMOBFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "FAULTMOB"));
RigFemScalarResultFrames * PCRITFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "PCRIT"));
frameCountProgress.incrementProgress();
const RigFemPart * femPart = m_femParts->part(partIndex);
const std::vector<cvf::Vec3f>& nodeCoordinates = femPart->nodes().coordinates;
float tanFricAng = tan(m_frictionAngleRad);
float cohPrTanFricAngle = (float)(m_cohesion/tanFricAng);
int frameCount = s11Frames->frameCount();
for(int fIdx = 0; fIdx < frameCount; ++fIdx)
{
@ -955,21 +973,30 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateSurfaceAlignedSt
const std::vector<float>& s23 = s23Frames->frameData(fIdx);
const std::vector<float>& s13 = s13Frames->frameData(fIdx);
std::vector<float>& sNorm = sNormFrames->frameData(fIdx);
std::vector<float>& sHori = sHoriFrames->frameData(fIdx);
std::vector<float>& sVert = sVertFrames->frameData(fIdx);
std::vector<float>& sNH = sNHFrames->frameData(fIdx);
std::vector<float>& sHV = sNVFrames->frameData(fIdx);
std::vector<float>& sNV = sHVFrames->frameData(fIdx);
std::vector<float>& SNDat = SNFrames->frameData(fIdx);
std::vector<float>& STHDat = STHFrames->frameData(fIdx);
std::vector<float>& STQVDat = STQVFrames->frameData(fIdx);
std::vector<float>& TNHDat = TNHFrames->frameData(fIdx);
std::vector<float>& TNQVDat = TNQVFrames->frameData(fIdx);
std::vector<float>& THQVDat = THQVFrames->frameData(fIdx);
std::vector<float>& TPDat = TPFrames->frameData(fIdx);
std::vector<float>& TincDat = TPincFrames->frameData(fIdx);
std::vector<float>& FAULTMOBDat = FAULTMOBFrames->frameData(fIdx);
std::vector<float>& PCRITDat = PCRITFrames->frameData(fIdx);
// HACK ! Todo : make it robust against other elements than Hex8
size_t valCount = s11.size() * 3; // Number of Elm Node Face results 24 = 4 * num faces = 3* numElmNodes
sNorm.resize(valCount);
sHori.resize(valCount);
sVert.resize(valCount);
sNH.resize(valCount);
sHV.resize(valCount);
sNV.resize(valCount);
SNDat.resize(valCount);
STHDat.resize(valCount);
STQVDat.resize(valCount);
TNHDat.resize(valCount);
TNQVDat.resize(valCount);
THQVDat.resize(valCount);
TPDat.resize(valCount);
TincDat.resize(valCount);
FAULTMOBDat.resize(valCount);
PCRITDat.resize(valCount);
int elementCount = femPart->elementCount();
for(int elmIdx = 0; elmIdx < elementCount; ++elmIdx)
@ -1017,12 +1044,31 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateSurfaceAlignedSt
caf::Ten3f xfTen = tensor.rotated(rotMx);
int elmNodFaceResIdx = elmNodFaceResIdxFaceStart + qIdx;
sHori[elmNodFaceResIdx]= xfTen[caf::Ten3f::SXX];
sVert[elmNodFaceResIdx]= xfTen[caf::Ten3f::SYY];
sNorm[elmNodFaceResIdx]= xfTen[caf::Ten3f::SZZ];
sNH[elmNodFaceResIdx]= xfTen[caf::Ten3f::SZX];
sHV[elmNodFaceResIdx]= xfTen[caf::Ten3f::SXY];
sNV[elmNodFaceResIdx]= xfTen[caf::Ten3f::SYZ];
float szx = xfTen[caf::Ten3f::SZX];
float syz = xfTen[caf::Ten3f::SYZ];
float szz = xfTen[caf::Ten3f::SZZ];
STHDat[elmNodFaceResIdx] = xfTen[caf::Ten3f::SXX];
STQVDat[elmNodFaceResIdx] = xfTen[caf::Ten3f::SYY];
SNDat[elmNodFaceResIdx] = xfTen[caf::Ten3f::SZZ];
TNHDat[elmNodFaceResIdx] = xfTen[caf::Ten3f::SZX];
TNQVDat[elmNodFaceResIdx] = xfTen[caf::Ten3f::SXY];
THQVDat[elmNodFaceResIdx] = xfTen[caf::Ten3f::SYZ];
float TP = sqrt(szx*szx + syz*syz);
TPDat[elmNodFaceResIdx] = TP;
if (TP > 1e-5)
{
TincDat[elmNodFaceResIdx] = cvf::Math::toDegrees( acos(syz/TP) );
}
else
{
TincDat[elmNodFaceResIdx] = std::numeric_limits<float>::infinity();
}
FAULTMOBDat[elmNodFaceResIdx] = TP/(tanFricAng * (szz + cohPrTanFricAngle));
PCRITDat[elmNodFaceResIdx] = szz - TP/tanFricAng;
}
}
}

View File

@ -47,6 +47,8 @@ public:
void setActiveFormationNames(RigFormationNames* activeFormationNames);
RigFormationNames* activeFormationNames();
void setCalculationParameters(double cohesion, double frictionAngleRad);
double parameterCohesion() const { return m_cohesion;}
double parameterFrictionAngleRad() const { return m_frictionAngleRad; }
std::map<std::string, std::vector<std::string> > scalarFieldAndComponentNames(RigFemResultPosEnum resPos);
std::vector<std::string> stepNames();

View File

@ -20,13 +20,64 @@
#include "RigFemPartResultsCollection.h"
#include "cvfGeometryTools.h"
#include "RivHexGridIntersectionTools.h"
#include "cvfAssert.h"
float RiuGeoMechXfTensorResultAccessor::SN (const caf::Ten3f& t) const { return t[caf::Ten3f::SZZ];}
float RiuGeoMechXfTensorResultAccessor::STH (const caf::Ten3f& t) const { return t[caf::Ten3f::SXX]; }
float RiuGeoMechXfTensorResultAccessor::STQV (const caf::Ten3f& t) const { return t[caf::Ten3f::SYY]; }
float RiuGeoMechXfTensorResultAccessor::TNH (const caf::Ten3f& t) const { return t[caf::Ten3f::SZX]; }
float RiuGeoMechXfTensorResultAccessor::TNQV (const caf::Ten3f& t) const { return t[caf::Ten3f::SYZ]; }
float RiuGeoMechXfTensorResultAccessor::THQV (const caf::Ten3f& t) const { return t[caf::Ten3f::SXY]; }
float RiuGeoMechXfTensorResultAccessor::TP (const caf::Ten3f& t) const
{
float szx = t[caf::Ten3f::SZX];
float szy = t[caf::Ten3f::SYZ];
float tp = sqrt(szx*szx + szy*szy);
return tp;
}
float RiuGeoMechXfTensorResultAccessor::TPinc (const caf::Ten3f& t) const
{
float szy = t[caf::Ten3f::SYZ];
float tp = TP(t);
if ( tp > 1e-5 )
{
return cvf::Math::toDegrees(acos(szy/tp));
}
else
{
return std::numeric_limits<float>::infinity();
}
}
float RiuGeoMechXfTensorResultAccessor::FAULTMOB(const caf::Ten3f& t) const
{
float szz = t[caf::Ten3f::SZZ];
float tp = TP(t);
return tp/(m_tanFricAng * (szz + m_cohPrTanFricAngle));
}
float RiuGeoMechXfTensorResultAccessor::PCRIT (const caf::Ten3f& t) const
{
float szz = t[caf::Ten3f::SZZ];
float tp = TP(t);
return szz - tp/m_tanFricAng;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RiuGeoMechXfTensorResultAccessor::RiuGeoMechXfTensorResultAccessor(RigFemPartResultsCollection * femResCollection, const RigFemResultAddress& resVarAddress, int timeStepIdx)
RiuGeoMechXfTensorResultAccessor::RiuGeoMechXfTensorResultAccessor(RigFemPartResultsCollection * femResCollection,
const RigFemResultAddress& resVarAddress,
int timeStepIdx)
{
RigFemResultAddress tensComp = resVarAddress;
tensComp.resultPosType = RIG_ELEMENT_NODAL;
@ -44,14 +95,19 @@ RiuGeoMechXfTensorResultAccessor::RiuGeoMechXfTensorResultAccessor(RigFemPartRes
tensComp.componentName = "S13";
tens13 = &femResCollection->resultValues(tensComp, 0, timeStepIdx);
resultComponent = caf::Ten3f::SZZ;
if ( resVarAddress.componentName == "SN" ) { m_tensorOperation = &RiuGeoMechXfTensorResultAccessor::SN ;}
if ( resVarAddress.componentName == "STH" ) { m_tensorOperation = &RiuGeoMechXfTensorResultAccessor::STH ;}
if ( resVarAddress.componentName == "STQV" ) { m_tensorOperation = &RiuGeoMechXfTensorResultAccessor::STQV ;}
if ( resVarAddress.componentName == "TNH" ) { m_tensorOperation = &RiuGeoMechXfTensorResultAccessor::TNH ;}
if ( resVarAddress.componentName == "TNQV" ) { m_tensorOperation = &RiuGeoMechXfTensorResultAccessor::TNQV ;}
if ( resVarAddress.componentName == "THQV" ) { m_tensorOperation = &RiuGeoMechXfTensorResultAccessor::THQV ;}
if ( resVarAddress.componentName == "TP" ) { m_tensorOperation = &RiuGeoMechXfTensorResultAccessor::TP ;}
if ( resVarAddress.componentName == "TPinc" ) { m_tensorOperation = &RiuGeoMechXfTensorResultAccessor::TPinc ;}
if ( resVarAddress.componentName == "FAULTMOB" ) { m_tensorOperation = &RiuGeoMechXfTensorResultAccessor::FAULTMOB;}
if ( resVarAddress.componentName == "PCRIT" ) { m_tensorOperation = &RiuGeoMechXfTensorResultAccessor::PCRIT ;}
if ( resVarAddress.componentName == "SN" ) resultComponent = caf::Ten3f::SZZ;
if ( resVarAddress.componentName == "STH" ) resultComponent = caf::Ten3f::SXX;
if ( resVarAddress.componentName == "STQV" ) resultComponent = caf::Ten3f::SYY;
if ( resVarAddress.componentName == "TNH" ) resultComponent = caf::Ten3f::SZX;
if ( resVarAddress.componentName == "TNQV" ) resultComponent = caf::Ten3f::SYZ;
if ( resVarAddress.componentName == "THQV" ) resultComponent = caf::Ten3f::SXY;
m_tanFricAng = tan((float)femResCollection->parameterFrictionAngleRad());
m_cohPrTanFricAngle = (float)((float)femResCollection->parameterCohesion()/m_tanFricAng);
}
//--------------------------------------------------------------------------------------------------
@ -104,7 +160,7 @@ void RiuGeoMechXfTensorResultAccessor::calculateInterpolatedValue(const cvf::Vec
ipT12, ipT23, ipT13);
caf::Ten3f xfTen = tensor.rotated(triangleXf);
returnValues[triangleVxIdx] = xfTen[resultComponent];
returnValues[triangleVxIdx] = m_tensorOperation(*this, xfTen);
}
}
}
@ -140,8 +196,7 @@ float RiuGeoMechXfTensorResultAccessor::calculateElmNodeValue(const std::array<c
ipT12, ipT23, ipT13);
caf::Ten3f xfTen = tensor.rotated(triangleXf);
float scalarValue = xfTen[resultComponent];
float scalarValue = m_tensorOperation(*this, xfTen);
return scalarValue;
}
}

View File

@ -21,6 +21,7 @@
#include <array>
#include "cafTensor3.h"
#include <functional>
class RigFemPartResultsCollection;
class RigFemResultAddress;
@ -29,7 +30,9 @@ class RivIntersectionVertexWeights;
class RiuGeoMechXfTensorResultAccessor
{
public:
RiuGeoMechXfTensorResultAccessor(RigFemPartResultsCollection * femResCollection, const RigFemResultAddress& resVarAddress, int timeStepIdx);
RiuGeoMechXfTensorResultAccessor(RigFemPartResultsCollection * femResCollection,
const RigFemResultAddress& resVarAddress,
int timeStepIdx);
void calculateInterpolatedValue(const cvf::Vec3f triangle[3], const RivIntersectionVertexWeights vertexWeights[3], float returnValues[3]);
@ -45,7 +48,21 @@ private:
const std::vector<float>* tens23;
const std::vector<float>* tens13;
caf::Ten3f::TensorComponentEnum resultComponent;
float m_tanFricAng ;
float m_cohPrTanFricAngle;
std::function<float (const RiuGeoMechXfTensorResultAccessor&, const caf::Ten3f&) > m_tensorOperation;
float SN (const caf::Ten3f& t) const;
float STH (const caf::Ten3f& t) const;
float STQV (const caf::Ten3f& t) const;
float TNH (const caf::Ten3f& t) const;
float TNQV (const caf::Ten3f& t) const;
float THQV (const caf::Ten3f& t) const;
float TP (const caf::Ten3f& t) const;
float TPinc (const caf::Ten3f& t) const;
float FAULTMOB(const caf::Ten3f& t) const;
float PCRIT (const caf::Ten3f& t) const;
};