Geom Range WIP

Store face in neighbor
Oposite face lookup
Start of finding grid corner element

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPart.cpp

@@ -154,11 +154,11 @@ void RigFemPart::calculateElmNeighbors()
         const int* elmNodes = this->connectivities(eIdx);
 
         int faceCount = RigFemTypes::elmentFaceCount(elmType);
-
+        int neighborCount = 0;
         for (int faceIdx = 0; faceIdx < faceCount; ++faceIdx)
         {
-            m_elmNeighbors[eIdx].idxToNeighborElmPrFace[faceIdx] = -1;
-
+            m_elmNeighbors[eIdx].indicesToNeighborElms[faceIdx] = -1;
+            m_elmNeighbors[eIdx].faceInNeighborElm[faceIdx] = -1;
             int faceNodeCount = 0;
             const int* localFaceIndices = RigFemTypes::localElmNodeIndicesForFace(elmType, faceIdx, &faceNodeCount);
 
@@ -222,15 +222,59 @@ void RigFemPart::calculateElmNeighbors()
                     // Compare faces
                     if (fComp.isSameButOposite(nbcElmNodes, nbcLocalFaceIndices, nbcFaceNodeCount))
                     {
-                        m_elmNeighbors[eIdx].idxToNeighborElmPrFace[faceIdx] = nbcElmIdx;
+                        m_elmNeighbors[eIdx].indicesToNeighborElms[faceIdx] = nbcElmIdx;
+                        m_elmNeighbors[eIdx].faceInNeighborElm[faceIdx] = nbcFaceIdx;
                         isNeighborFound = true;
+                       
                         break;
                     }
                 }
 
-                if (isNeighborFound) break;
+                if (isNeighborFound)
+                {
+                    ++neighborCount;
+                    break;
+                }
             }
         }
+
+        if ((faceCount - neighborCount) >= 3)
+        {
+            m_possibleGridCornerElements.push_back(eIdx);
+        }
     }
 }
 
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+cvf::Vec3f RigFemPart::faceNormal(int elmIdx, int faceIdx)
+{
+    const std::vector<cvf::Vec3f>& nodeCoordinates = this->nodes().coordinates;
+
+    RigElementType eType = this->elementType(elmIdx);
+    const int* elmNodeIndices = this->connectivities(elmIdx);
+
+    int faceNodeCount = 0;
+    const int*  localElmNodeIndicesForFace = RigFemTypes::localElmNodeIndicesForFace(eType, faceIdx, &faceNodeCount);
+    
+    if (faceNodeCount == 4)
+    {
+        const 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::Vec3f normal = (*(quadVxs[2]) -  *(quadVxs[0])) ^ (*(quadVxs[3]) - *(quadVxs[1]));
+        return normal;
+    }
+    else if (faceNodeCount != 4)
+    {
+        CVF_ASSERT(false);
+    }
+
+    return cvf::Vec3f::ZERO;
+}
+

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPart.h

@@ -67,7 +67,11 @@ public:
     
     void                        assertElmNeighborsIsCalculated();
     int                         elementNeighbor(int elementIndex, int faceIndex) const
-                                { return m_elmNeighbors[elementIndex].idxToNeighborElmPrFace[faceIndex]; }
+                                { return m_elmNeighbors[elementIndex].indicesToNeighborElms[faceIndex]; }
+    const std::vector<int>&     possibleGridCornerElements() const { return m_possibleGridCornerElements; }
+
+    cvf::Vec3f                  faceNormal(int elmentIndex, int faceIndex);
+
     const RigFemPartGrid*       structGrid();   
 
 private:
@@ -86,6 +90,8 @@ private:
     std::vector<std::vector<size_t> > m_nodeToElmRefs; // Needs a more memory friendly structure
   
     void calculateElmNeighbors();
-    struct Neighbors { int idxToNeighborElmPrFace[6]; };
+    struct Neighbors { int indicesToNeighborElms[6]; char faceInNeighborElm[6];};
     std::vector<  Neighbors > m_elmNeighbors;
+    std::vector<int> m_possibleGridCornerElements;
+
 };

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartGrid.cpp

@@ -40,6 +40,33 @@ RigFemPartGrid::~RigFemPartGrid()
 
 }
 
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+
+cvf::StructGridInterface::FaceType RigFemPartGrid::findGridFace(cvf::Vec3d faceNormal )
+{
+    FaceType bestFace = cvf::StructGridInterface::POS_I;
+
+    double maxComponent = fabs(faceNormal[0]);
+    bestFace = (faceNormal[0] < 0) ? cvf::StructGridInterface::NEG_I: cvf::StructGridInterface::POS_I;
+
+    double absComp = fabs(faceNormal[1]);
+    if ( absComp > maxComponent)
+    {
+        maxComponent = absComp;
+        bestFace = (faceNormal[1] < 0) ? cvf::StructGridInterface::NEG_J: cvf::StructGridInterface::POS_J;
+    }
+
+    absComp = fabs(faceNormal[2]);
+    if ( absComp > maxComponent)
+    {
+        bestFace = (faceNormal[2] < 0) ? cvf::StructGridInterface::NEG_K: cvf::StructGridInterface::POS_K;
+    }
+
+    return bestFace;
+}
+
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
@@ -62,10 +89,30 @@ void RigFemPartGrid::generateStructGridData()
     // 6. If IJK to elm idx is needed, allocate "grid" with maxI,maxJ,maxZ values
     //    Loop over elms, assign elmIdx to IJK address in grid
     
+    const std::vector<int>& possibleGridCorners = m_femPart->possibleGridCornerElements();
+    size_t possibleCornerCount = possibleGridCorners.size();
+    const std::vector<cvf::Vec3f>& nodeCoordinates = m_femPart->nodes().coordinates;
+
+    // Find corner cell closest to origo
+    size_t gridCornerClosestToOrigo = cvf::UNDEFINED_SIZE_T;
+    double minDistance =  HUGE_VAL;
+    for (size_t pcIdx = 0; pcIdx < possibleCornerCount; ++pcIdx)
+    {
+        int elmIdx =  possibleGridCorners[pcIdx];
+
+        const int* elmNodeIndices =  m_femPart->connectivities(elmIdx);
+        cvf::Vec3f firstNodePos = nodeCoordinates[elmNodeIndices[0]];
+        float distSq = firstNodePos.lengthSquared();
+        if (distSq < minDistance)
+        {
+            minDistance = distSq;
+            gridCornerClosestToOrigo = pcIdx;
+        }
+    }
+
 
 }
 
-
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------

ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartGrid.h

@@ -50,7 +50,7 @@ public:
  
  private:
     void generateStructGridData();
-
+    static FaceType findGridFace(cvf::Vec3d faceNormal);
 
     RigFemPart* m_femPart;
 };

ApplicationCode/GeoMech/GeoMechDataModel/RigFemTypes.cpp

@@ -77,3 +77,23 @@ const int* RigFemTypes::localElmNodeIndicesForFace(RigElementType elmType, int f
 
     return CAX4_Faces;
 }
+
+int RigFemTypes::oppositeFace(RigElementType elmType, int faceIdx)
+{
+    static const int HEX8_OppositeFaces[6]={ 1, 0, 3, 2, 5, 4};
+ 
+    switch (elmType)
+    {
+        case HEX8:
+            return HEX8_OppositeFaces[faceIdx];
+            break;
+        case CAX4:
+            return faceIdx;
+            break;
+        default:
+            assert(false); // Element type not supported
+            break;
+    }
+
+    return faceIdx;
+}

ApplicationCode/GeoMech/GeoMechDataModel/RigFemTypes.h

@@ -33,5 +33,5 @@ public:
     static int elmentNodeCount(RigElementType elmType);
     static int elmentFaceCount(RigElementType elmType);
     static const int* localElmNodeIndicesForFace(RigElementType elmType, int faceIdx, int* faceNodeCount);
-
+    static int opositeFace(RigElementType elmType, int faceIdx);
 };