#942 Added plane Azi and Inc angles as result values.

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartResultsCollection.cpp

@@ -339,6 +339,8 @@ std::map<std::string, std::vector<std::string> > 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<cvf::Vec3f>& 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<float>& Pazi = PaziFrames->frameData(fIdx);
+        std::vector<float>& 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,6 +1365,9 @@ RigFemScalarResultFrames* RigFemPartResultsCollection::calculateDerivedResult(in
 
     if(resVarAddr.resultPosType == RIG_ELEMENT_NODAL_FACE )
     {
+        if (resVarAddr.componentName == "Pazi" || resVarAddr.componentName == "Pinc" )
+            return calculateSurfaceAngles(partIndex, resVarAddr);
+        else
             return calculateSurfaceAlignedStress(partIndex, resVarAddr);
     }
 

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 inc() const { return incAziR[0];}
-    float azi() { return incAziR[1];}
+    float azi() const { return incAziR[1];}
-    float r()   { return incAziR[2];}
+    float r()   const { return incAziR[2];}
 
 private:
     std::array<float, 3> incAziR;

ApplicationCode/ModelVisualization/Intersections/RivIntersectionBoxPartMgr.cpp

@@ -184,6 +184,15 @@ void RivIntersectionBoxPartMgr::updateCellResultColor(size_t timeStepIndex)
             // Special direction sensitive result calculation
             const cvf::Vec3fArray* triangelVxes = m_intersectionBoxGenerator->triangleVxes();
 
+            if ( resVarAddress.componentName == "Pazi" || resVarAddress.componentName == "Pinc" )
+            {
+                RivIntersectionPartMgr::calculatePlaneAngleTextureCoords(m_intersectionBoxFacesTextureCoords.p(),
+                                                                         triangelVxes,
+                                                                         resVarAddress,
+                                                                         mapper);
+            }
+            else
+            {
                 RivIntersectionPartMgr::calculateGeoMechTensorXfTextureCoords(m_intersectionBoxFacesTextureCoords.p(),
                                                                           triangelVxes,
                                                                           vertexWeights,
@@ -192,6 +201,7 @@ void RivIntersectionBoxPartMgr::updateCellResultColor(size_t timeStepIndex)
                                                                           (int)timeStepIndex,
                                                                           mapper);
             }
+        }
 
         RivScalarMapperUtils::applyTextureResultsToPart(m_intersectionBoxFaces.p(), 
                                                         m_intersectionBoxFacesTextureCoords.p(), 

ApplicationCode/ModelVisualization/Intersections/RivIntersectionPartMgr.cpp

@@ -188,6 +188,15 @@ void RivIntersectionPartMgr::updateCellResultColor(size_t timeStepIndex)
             // Special direction sensitive result calculation
             const cvf::Vec3fArray* triangelVxes = m_crossSectionGenerator->triangleVxes();
 
+            if ( resVarAddress.componentName == "Pazi" || resVarAddress.componentName == "Pinc" )
+            {
+                RivIntersectionPartMgr::calculatePlaneAngleTextureCoords(m_crossSectionFacesTextureCoords.p(),
+                                                                         triangelVxes, 
+                                                                         resVarAddress,
+                                                                         mapper);
+            }
+            else
+            {
                 RivIntersectionPartMgr::calculateGeoMechTensorXfTextureCoords(m_crossSectionFacesTextureCoords.p(),
                                                                               triangelVxes,
                                                                               vertexWeights,
@@ -196,6 +205,7 @@ void RivIntersectionPartMgr::updateCellResultColor(size_t timeStepIndex)
                                                                               (int)timeStepIndex,
                                                                               mapper);
             }
+        }
 
         RivScalarMapperUtils::applyTextureResultsToPart(m_crossSectionFaces.p(), 
                                                         m_crossSectionFacesTextureCoords.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<int>(triangelVertices->size());
+    int triCount = vxCount/3;
+
+    std::function<float (const OffshoreSphericalCoords& )> 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<float>::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

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<RivIntersectionHexGridInterface> createHexGridInterface();
 private: