From d2d9d43ca14d4dde8e81cd03ce3b99d7aa14523b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jacob=20St=C3=B8ren?= Date: Tue, 1 Nov 2016 12:04:56 +0100 Subject: [PATCH] #942 Added plane Azi and Inc angles as result values. --- .../RigFemPartResultsCollection.cpp | 83 ++++++++++++++++++- .../RigFemPartResultsCollection.h | 7 +- .../RivIntersectionBoxPartMgr.cpp | 12 ++- .../Intersections/RivIntersectionPartMgr.cpp | 69 +++++++++++++-- .../Intersections/RivIntersectionPartMgr.h | 5 ++ 5 files changed, 162 insertions(+), 14 deletions(-) diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp index 64231f9f7c..14ae2c3b55 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp @@ -339,6 +339,8 @@ std::map > RigFemPartResultsCollection::sc } else if(resPos == RIG_ELEMENT_NODAL_FACE) { + fieldCompNames["Plane"].push_back("Pinc"); + fieldCompNames["Plane"].push_back("Pazi"); fieldCompNames["SE"].push_back("SN"); fieldCompNames["SE"].push_back("STH"); @@ -631,7 +633,6 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateSFI(int partInde return dstDataFrames; } -#define M_PI_4 0.785398163397448309616 // pi/4 //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -670,7 +671,8 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateDSM(int partInde { float se1 = se1Data[vIdx]; float se3 = se3Data[vIdx]; - float rho = 2.0f * atan( sqrt(( se1 + cohPrTanFricAngle)/(se3 + cohPrTanFricAngle)) - M_PI_4); + float pi_4 = 0.785398163397448309616f; + float rho = 2.0f * atan( sqrt(( se1 + cohPrTanFricAngle)/(se3 + cohPrTanFricAngle)) - pi_4); { dstFrameData[vIdx] = tan(rho)/tanFricAng; @@ -1081,6 +1083,78 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateSurfaceAlignedSt return requestedSurfStress; } +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +RigFemScalarResultFrames* RigFemPartResultsCollection::calculateSurfaceAngles(int partIndex, const RigFemResultAddress& resVarAddr) +{ + CVF_ASSERT(resVarAddr.componentName == "Pazi" || resVarAddr.componentName == "Pinc"); + + caf::ProgressInfo frameCountProgress(this->frameCount() * 1, ""); + frameCountProgress.setProgressDescription("Calculating " + QString::fromStdString(resVarAddr.fieldName + ": " + resVarAddr.componentName)); + + RigFemScalarResultFrames * PaziFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "Pazi")); + RigFemScalarResultFrames * PincFrames = m_femPartResults[partIndex]->createScalarResult(RigFemResultAddress(resVarAddr.resultPosType, resVarAddr.fieldName, "Pinc")); + + const RigFemPart * femPart = m_femParts->part(partIndex); + const std::vector& nodeCoordinates = femPart->nodes().coordinates; + int frameCount = this->frameCount(); + + // HACK ! Todo : make it robust against other elements than Hex8 + size_t valCount = femPart->elementCount() * 24; // Number of Elm Node Face results 24 = 4 * num faces = 3* numElmNodes + + for ( int fIdx = 0; fIdx < frameCount; ++fIdx ) + { + std::vector& Pazi = PaziFrames->frameData(fIdx); + std::vector& Pinc = PincFrames->frameData(fIdx); + + Pazi.resize(valCount); + Pinc.resize(valCount); + + int elementCount = femPart->elementCount(); + for ( int elmIdx = 0; elmIdx < elementCount; ++elmIdx ) + { + RigElementType elmType = femPart->elementType(elmIdx); + int faceCount = RigFemTypes::elmentFaceCount(elmType); + const int* elmNodeIndices = femPart->connectivities(elmIdx); + + int elmNodFaceResIdxElmStart = elmIdx * 24; // HACK should get from part + + for ( int lfIdx = 0; lfIdx < faceCount; ++lfIdx ) + { + int faceNodeCount = 0; + const int* localElmNodeIndicesForFace = RigFemTypes::localElmNodeIndicesForFace(elmType, lfIdx, &faceNodeCount); + if ( faceNodeCount == 4 ) + { + int elmNodFaceResIdxFaceStart = elmNodFaceResIdxElmStart + lfIdx*4; // HACK + cvf::Vec3f quadVxs[4]; + + quadVxs[0] = (nodeCoordinates[elmNodeIndices[localElmNodeIndicesForFace[0]]]); + quadVxs[1] = (nodeCoordinates[elmNodeIndices[localElmNodeIndicesForFace[1]]]); + quadVxs[2] = (nodeCoordinates[elmNodeIndices[localElmNodeIndicesForFace[2]]]); + quadVxs[3] = (nodeCoordinates[elmNodeIndices[localElmNodeIndicesForFace[3]]]); + + cvf::Mat3f rotMx = cvf::GeometryTools::computePlaneHorizontalRotationMx(quadVxs[2] -quadVxs[0], quadVxs[3] - quadVxs[1]); + OffshoreSphericalCoords sphCoord(cvf::Vec3f(rotMx.rowCol(0,2), rotMx.rowCol(1,2), rotMx.rowCol(2,2))); // Use Ez from the matrix as plane normal + + for ( int qIdx = 0; qIdx < 4; ++qIdx ) + { + int elmNodFaceResIdx = elmNodFaceResIdxFaceStart + qIdx; + Pazi[elmNodFaceResIdx] = cvf::Math::toDegrees( sphCoord.azi() ); + Pinc[elmNodFaceResIdx] = cvf::Math::toDegrees( sphCoord.inc() ); + } + } + } + } + + frameCountProgress.incrementProgress(); + } + + RigFemScalarResultFrames* requestedPlaneAngle = this->findOrLoadScalarResult(partIndex, resVarAddr); + return requestedPlaneAngle; + +} + //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- @@ -1291,7 +1365,10 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateDerivedResult(in if(resVarAddr.resultPosType == RIG_ELEMENT_NODAL_FACE ) { - return calculateSurfaceAlignedStress(partIndex, resVarAddr); + if (resVarAddr.componentName == "Pazi" || resVarAddr.componentName == "Pinc" ) + return calculateSurfaceAngles(partIndex, resVarAddr); + else + return calculateSurfaceAlignedStress(partIndex, resVarAddr); } if (resVarAddr.fieldName == "SE" && resVarAddr.componentName == "SFI") diff --git a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h index 55803c5fde..1898abb026 100644 --- a/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h +++ b/ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.h @@ -98,6 +98,7 @@ private: RigFemScalarResultFrames* calculateVolumetricStrain(int partIndex, const RigFemResultAddress& resVarAddr); RigFemScalarResultFrames* calculateDeviatoricStrain(int partIndex, const RigFemResultAddress& resVarAddr); RigFemScalarResultFrames* calculateSurfaceAlignedStress(int partIndex, const RigFemResultAddress& resVarAddr); + RigFemScalarResultFrames* calculateSurfaceAngles(int partIndex, const RigFemResultAddress& resVarAddr); RigFemScalarResultFrames* calculatePrincipalStressValues(int partIndex, const RigFemResultAddress &resVarAddr); RigFemScalarResultFrames* calculatePrincipalStrainValues(int partIndex, const RigFemResultAddress &resVarAddr); @@ -140,9 +141,9 @@ public: } - float inc() { return incAziR[0];} - float azi() { return incAziR[1];} - float r() { return incAziR[2];} + float inc() const { return incAziR[0];} + float azi() const { return incAziR[1];} + float r() const { return incAziR[2];} private: std::array incAziR; diff --git a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionBoxPartMgr.cpp b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionBoxPartMgr.cpp index 147871d2e9..d5d8ca1755 100644 --- a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionBoxPartMgr.cpp +++ b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionBoxPartMgr.cpp @@ -184,13 +184,23 @@ void RivIntersectionBoxPartMgr::updateCellResultColor(size_t timeStepIndex) // Special direction sensitive result calculation const cvf::Vec3fArray* triangelVxes = m_intersectionBoxGenerator->triangleVxes(); - RivIntersectionPartMgr::calculateGeoMechTensorXfTextureCoords(m_intersectionBoxFacesTextureCoords.p(), + if ( resVarAddress.componentName == "Pazi" || resVarAddress.componentName == "Pinc" ) + { + RivIntersectionPartMgr::calculatePlaneAngleTextureCoords(m_intersectionBoxFacesTextureCoords.p(), + triangelVxes, + resVarAddress, + mapper); + } + else + { + RivIntersectionPartMgr::calculateGeoMechTensorXfTextureCoords(m_intersectionBoxFacesTextureCoords.p(), triangelVxes, vertexWeights, caseData, resVarAddress, (int)timeStepIndex, mapper); + } } RivScalarMapperUtils::applyTextureResultsToPart(m_intersectionBoxFaces.p(), diff --git a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionPartMgr.cpp b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionPartMgr.cpp index ca80711b39..64bf8d9f53 100644 --- a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionPartMgr.cpp +++ b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionPartMgr.cpp @@ -188,13 +188,23 @@ void RivIntersectionPartMgr::updateCellResultColor(size_t timeStepIndex) // Special direction sensitive result calculation const cvf::Vec3fArray* triangelVxes = m_crossSectionGenerator->triangleVxes(); - RivIntersectionPartMgr::calculateGeoMechTensorXfTextureCoords(m_crossSectionFacesTextureCoords.p(), - triangelVxes, - vertexWeights, - caseData, - resVarAddress, - (int)timeStepIndex, - mapper); + if ( resVarAddress.componentName == "Pazi" || resVarAddress.componentName == "Pinc" ) + { + RivIntersectionPartMgr::calculatePlaneAngleTextureCoords(m_crossSectionFacesTextureCoords.p(), + triangelVxes, + resVarAddress, + mapper); + } + else + { + RivIntersectionPartMgr::calculateGeoMechTensorXfTextureCoords(m_crossSectionFacesTextureCoords.p(), + triangelVxes, + vertexWeights, + caseData, + resVarAddress, + (int)timeStepIndex, + mapper); + } } RivScalarMapperUtils::applyTextureResultsToPart(m_crossSectionFaces.p(), @@ -287,6 +297,51 @@ void RivIntersectionPartMgr::calculateGeoMechTensorXfTextureCoords(cvf::Vec2fArr } + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +void RivIntersectionPartMgr::calculatePlaneAngleTextureCoords(cvf::Vec2fArray* textureCoords, + const cvf::Vec3fArray* triangelVertices, + const RigFemResultAddress& resVarAddress, + const cvf::ScalarMapper* mapper) +{ + + textureCoords->resize(triangelVertices->size()); + cvf::Vec2f* rawPtr = textureCoords->ptr(); + int vxCount = static_cast(triangelVertices->size()); + int triCount = vxCount/3; + + std::function operation; + if (resVarAddress.componentName == "Pazi") + { + operation = [](const OffshoreSphericalCoords& sphCoord) { return sphCoord.azi();}; + } + else if ( resVarAddress.componentName == "Pinc" ) + { + operation = [](const OffshoreSphericalCoords& sphCoord) { return sphCoord.inc();}; + } + + #pragma omp parallel for schedule(dynamic) + for ( int triangleIdx = 0; triangleIdx < triCount; ++triangleIdx ) + { + int triangleVxStartIdx = triangleIdx*3; + + const cvf::Vec3f* triangle = &((*triangelVertices)[triangleVxStartIdx]); + cvf::Mat3f rotMx = cvf::GeometryTools::computePlaneHorizontalRotationMx(triangle[1] - triangle[0], triangle[2] - triangle[0]); + + OffshoreSphericalCoords sphCoord(cvf::Vec3f(rotMx.rowCol(0, 2), rotMx.rowCol(1, 2), rotMx.rowCol(2, 2))); // Use Ez from the matrix as plane normal + + float angle = cvf::Math::toDegrees( operation(sphCoord)); + cvf::Vec2f texCoord = (angle != std::numeric_limits::infinity()) ? mapper->mapToTextureCoord(angle) : cvf::Vec2f(0.0f, 1.0f); + rawPtr[triangleVxStartIdx + 0] = texCoord; + rawPtr[triangleVxStartIdx + 1] = texCoord; + rawPtr[triangleVxStartIdx + 2] = texCoord; + } + +} + + //-------------------------------------------------------------------------------------------------- /// Calculates the texture coordinates in a "nearly" one dimensional texture. /// Undefined values are coded with a y-texturecoordinate value of 1.0 instead of the normal 0.5 diff --git a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionPartMgr.h b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionPartMgr.h index 0da5083e0e..761ce01631 100644 --- a/ApplicationCode/ModelVisualization/Intersections/RivIntersectionPartMgr.h +++ b/ApplicationCode/ModelVisualization/Intersections/RivIntersectionPartMgr.h @@ -89,6 +89,11 @@ public: int timeStepIdx, const cvf::ScalarMapper* mapper); + static void calculatePlaneAngleTextureCoords(cvf::Vec2fArray* textureCoords, + const cvf::Vec3fArray* triangelVertices, + const RigFemResultAddress& resVarAddress, + const cvf::ScalarMapper* mapper); + cvf::ref createHexGridInterface(); private: