OilfieldToolsNet/ResInsight
4
0
Fork 0
mirror of https://github.com/OPM/ResInsight.git synced 2025-02-25 18:55:39 -06:00
Code Issues Packages Projects Releases Wiki Activity

Geomech Range filters working

Browse Source 
The algorithm for the IJK assignment now works on the simple test
example.
This commit is contained in:
Jacob Støren
2015-05-29 07:53:03 +02:00
parent 8a6e1ae65a
commit 51fd1b4de2
7 changed files with 305 additions and 81 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
ApplicationCode/GeoMech/GeoMechDataModel/RigFemPart.cpp
View File

@@ -97,7 +97,7 @@ void RigFemPart::calculateNodeToElmRefs()
{
m_nodeToElmRefs.resize(nodes().nodeIds.size());
for (size_t eIdx = 0; eIdx < m_elementId.size(); ++eIdx)
for (int eIdx = 0; eIdx < static_cast<int>(m_elementId.size()); ++eIdx)
{
int elmNodeCount = RigFemTypes::elmentNodeCount(elementType(eIdx));
const int* elmNodes = connectivities(eIdx);
@@ -111,7 +111,7 @@ void RigFemPart::calculateNodeToElmRefs()
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
const size_t* RigFemPart::elementsUsingNode(int nodeIndex)
const int* RigFemPart::elementsUsingNode(int nodeIndex)
{
return &(m_nodeToElmRefs[nodeIndex][0]);
}
@@ -167,7 +167,7 @@ void RigFemPart::calculateElmNeighbors()
{
int firstNodeIdxOfFace = elmNodes[localFaceIndices[0]];
int candidateCount1 = this->numElementsUsingNode(firstNodeIdxOfFace);
const size_t* candidates1 = this->elementsUsingNode(firstNodeIdxOfFace);
const int* candidates1 = this->elementsUsingNode(firstNodeIdxOfFace);
if (candidateCount1)
{
@@ -175,7 +175,7 @@ void RigFemPart::calculateElmNeighbors()
int thirdNodeIdxOfFace = elmNodes[localFaceIndices[3]];
int candidateCount2 = this->numElementsUsingNode(thirdNodeIdxOfFace);
const size_t* candidates2 = this->elementsUsingNode(thirdNodeIdxOfFace);
const int* candidates2 = this->elementsUsingNode(thirdNodeIdxOfFace);
// The candidates are sorted from smallest to largest, so we do a linear search to find the
// (two) common cells in the two arrays, and leaving this element out, we have one candidate left

7
ApplicationCode/GeoMech/GeoMechDataModel/RigFemPart.h
View File

@@ -62,12 +62,15 @@ public:
const RigFemPartNodes& nodes() const {return m_nodes;}
void assertNodeToElmIndicesIsCalculated();
const size_t* elementsUsingNode(int nodeIndex);
const int* elementsUsingNode(int nodeIndex);
int numElementsUsingNode(int nodeIndex);
void assertElmNeighborsIsCalculated();
int elementNeighbor(int elementIndex, int faceIndex) const
{ return m_elmNeighbors[elementIndex].indicesToNeighborElms[faceIndex]; }
int neighborFace(int elementIndex, int faceIndex) const
{ return m_elmNeighbors[elementIndex].faceInNeighborElm[faceIndex]; }
const std::vector<int>& possibleGridCornerElements() const { return m_possibleGridCornerElements; }
cvf::Vec3f faceNormal(int elmentIndex, int faceIndex);
@@ -87,7 +90,7 @@ private:
cvf::ref<RigFemPartGrid> m_structGrid;
void calculateNodeToElmRefs();
std::vector<std::vector<size_t> > m_nodeToElmRefs; // Needs a more memory friendly structure
std::vector<std::vector<int> > m_nodeToElmRefs; // Needs a more memory friendly structure
void calculateElmNeighbors();
struct Neighbors { int indicesToNeighborElms[6]; char faceInNeighborElm[6];};

327
ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartGrid.cpp
View File

@@ -43,83 +43,275 @@ RigFemPartGrid::~RigFemPartGrid()
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::StructGridInterface::FaceType RigFemPartGrid::findGridFace(cvf::Vec3d faceNormal )
void RigFemPartGrid::generateStructGridData()
{
FaceType bestFace = cvf::StructGridInterface::POS_I;
//[X] 1. Calculate neighbors for each element
//[X] record the ones with 3 or fewer neighbors as possible grid corners
//[X] 2. Loop over the possible corner cells,
//[X] find the one that corresponds to IJK = 000
//[X] by finding the one closest to origo
//[X] 4. Assign IJK = 000 to that element
//[X] Store IJK in elm idx array
//[X] 5. Loop along POS I surfaces increment I for each element and assign IJK
//[X] when at end, go to POS J neighbor, increment J, repeat above.
//[X] etc for POS Z
//[X] Find max IJK as you go,
//[ ] also assert that there are no NEG I/NEG J/NEG Z neighbors when starting on a new row
//[ ] (Need to find min, and offset IJK values if there exists such)
//[ ] 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
int gridCornerClosestToOrigo = findElmIdxOfGridCornerClosestToOrigo();
if (gridCornerClosestToOrigo == cvf::UNDEFINED_SIZE_T) return;
double maxComponent = fabs(faceNormal[0]);
bestFace = (faceNormal[0] < 0) ? cvf::StructGridInterface::NEG_I: cvf::StructGridInterface::POS_I;
// Find the IJK faces based on the corner cell
cvf::Vec3i ijkMainFaceIdx = cvf::Vec3i(-1,-1,-1);
double absComp = fabs(faceNormal[1]);
if ( absComp > maxComponent)
{
maxComponent = absComp;
bestFace = (faceNormal[1] < 0) ? cvf::StructGridInterface::NEG_J: cvf::StructGridInterface::POS_J;
RigElementType eType = m_femPart->elementType(gridCornerClosestToOrigo);
int faceCount = RigFemTypes::elmentFaceCount(eType);
std::vector<cvf::Vec3f> normals(faceCount);
for (int faceIdx = 0; faceIdx < faceCount; ++faceIdx)
{
normals[faceIdx] = m_femPart->faceNormal(gridCornerClosestToOrigo, faceIdx);
}
// Record three independent main direction vectors for the element, and what face they are created from
cvf::Vec3f mainElmDirections[3];
int mainElmDirOriginFaces[3];
if (eType == HEX8)
{
mainElmDirections[0] = normals[0] - normals[1]; // To get a better "average" direction vector
mainElmDirections[1] = normals[2] - normals[3];
mainElmDirections[2] = normals[4] - normals[5];
mainElmDirOriginFaces[0] = 0;
mainElmDirOriginFaces[1] = 2;
mainElmDirOriginFaces[2] = 4;
}
else
{
CVF_ASSERT(false);
}
// Match the element main directions with best XYZ match (IJK respectively)
// Find the max component of a mainElmDirection.
// Assign the index of that mainElmDirection to the mainElmDirectionIdxForIJK at the index of the max component.
int mainElmDirectionIdxForIJK[3] = { -1, -1, -1};
for (int dIdx = 0; dIdx < 3; ++dIdx)
{
double maxAbsComp = 0;
for (int cIdx = 2; cIdx >= 0 ; --cIdx)
{
float absComp = fabs(mainElmDirections[dIdx][cIdx]);
if (absComp > maxAbsComp)
{
maxAbsComp = absComp;
mainElmDirectionIdxForIJK[cIdx] = dIdx;
}
}
}
// make sure all the main directions are used
bool mainDirsUsed[3] = { false, false, false};
mainDirsUsed[mainElmDirectionIdxForIJK[0]] = true;
mainDirsUsed[mainElmDirectionIdxForIJK[1]] = true;
mainDirsUsed[mainElmDirectionIdxForIJK[2]] = true;
int unusedDir = -1;
if (!mainDirsUsed[0]) unusedDir = 0;
if (!mainDirsUsed[1]) unusedDir = 1;
if (!mainDirsUsed[2]) unusedDir = 2;
if (unusedDir >= 0)
{
if (mainElmDirectionIdxForIJK[0] == mainElmDirectionIdxForIJK[1]) mainElmDirectionIdxForIJK[0] = unusedDir;
else if (mainElmDirectionIdxForIJK[1] == mainElmDirectionIdxForIJK[2]) mainElmDirectionIdxForIJK[1] = unusedDir;
else if (mainElmDirectionIdxForIJK[2] == mainElmDirectionIdxForIJK[0]) mainElmDirectionIdxForIJK[2] = unusedDir;
}
// Assign the correct face based on the main direction
ijkMainFaceIdx[0] = (mainElmDirections[mainElmDirectionIdxForIJK[0]] * cvf::Vec3f::X_AXIS > 0) ? mainElmDirOriginFaces[mainElmDirectionIdxForIJK[0]]: RigFemTypes::oppositeFace(eType, mainElmDirOriginFaces[mainElmDirectionIdxForIJK[0]]);
ijkMainFaceIdx[1] = (mainElmDirections[mainElmDirectionIdxForIJK[1]] * cvf::Vec3f::Y_AXIS > 0) ? mainElmDirOriginFaces[mainElmDirectionIdxForIJK[1]]: RigFemTypes::oppositeFace(eType, mainElmDirOriginFaces[mainElmDirectionIdxForIJK[1]]);
ijkMainFaceIdx[2] = (mainElmDirections[mainElmDirectionIdxForIJK[2]] * -cvf::Vec3f::Z_AXIS > 0) ? mainElmDirOriginFaces[mainElmDirectionIdxForIJK[2]]: RigFemTypes::oppositeFace(eType, mainElmDirOriginFaces[mainElmDirectionIdxForIJK[2]]);
}
absComp = fabs(faceNormal[2]);
if ( absComp > maxComponent)
// assign ijk to cells
{
bestFace = (faceNormal[2] < 0) ? cvf::StructGridInterface::NEG_K: cvf::StructGridInterface::POS_K;
m_ijkPrElement.resize(m_femPart->elementCount(), cvf::Vec3i(-1,-1,-1));
int posIFaceIdx = ijkMainFaceIdx[0];
int posJFaceIdx = ijkMainFaceIdx[1];
int posKFaceIdx = ijkMainFaceIdx[2];
m_elmentIJKCounts = cvf::Vec3st(0, 0, 0);
int elmIdxInK = gridCornerClosestToOrigo;
cvf::Vec3f posKNormal = m_femPart->faceNormal(elmIdxInK, posKFaceIdx);
int kCoord = 0;
while (true)
{
int elmIdxInJ = elmIdxInK;
cvf::Vec3f startElmInKNormalJ = m_femPart->faceNormal(elmIdxInJ, posJFaceIdx);
cvf::Vec3f startElmInKNormalI = m_femPart->faceNormal(elmIdxInJ, posIFaceIdx);
int jCoord = 0;
while (true)
{
int elmIdxInI = elmIdxInJ;
cvf::Vec3f startElmInJNormalI = m_femPart->faceNormal(elmIdxInI, posIFaceIdx);
int iCoord = 0;
while (true)
{
// Assign ijk coordinate
m_ijkPrElement[elmIdxInI] = cvf::Vec3i(iCoord, jCoord, kCoord);
++iCoord;
// Find neighbor and exit if at end
int neighborElmIdx = m_femPart->elementNeighbor(elmIdxInI, posIFaceIdx);
if (neighborElmIdx == -1) break;
// Find the continuing face in the neighbor element (opposite of the neighbor face)
int neighborNegFaceIdx = m_femPart->neighborFace(elmIdxInI, posIFaceIdx);
RigElementType eType = m_femPart->elementType(neighborElmIdx);
posIFaceIdx = RigFemTypes::oppositeFace(eType, neighborNegFaceIdx);
// Step to neighbor
elmIdxInI = neighborElmIdx;
}
// Scoped to show that nothing bleeds further to K-loop
{
if (iCoord > m_elmentIJKCounts[0]) m_elmentIJKCounts[0] = iCoord;
++jCoord;
// Find neighbor and exit if at end
int neighborElmIdx = m_femPart->elementNeighbor(elmIdxInJ, posJFaceIdx);
if (neighborElmIdx == -1) break;
// Find the continuing face in the neighbor element (opposite of the neighbor face)
int neighborNegFaceIdx = m_femPart->neighborFace(elmIdxInJ, posJFaceIdx);
RigElementType eType = m_femPart->elementType(neighborElmIdx);
posJFaceIdx = RigFemTypes::oppositeFace(eType, neighborNegFaceIdx);
// Now where is posIFace of the new J cell ?
posIFaceIdx = perpendicularFaceInDirection(startElmInJNormalI, neighborNegFaceIdx, neighborElmIdx);
// Step to neighbor
elmIdxInJ = neighborElmIdx;
}
}
{
if (jCoord > m_elmentIJKCounts[1]) m_elmentIJKCounts[1] = jCoord;
++kCoord;
// Find neighbor and exit if at end
int neighborElmIdx = m_femPart->elementNeighbor(elmIdxInK, posKFaceIdx);
if (neighborElmIdx == -1) break;
// Find the continuing face in the neighbor element (opposite of the neighbor face)
int neighborNegFaceIdx = m_femPart->neighborFace(elmIdxInK, posKFaceIdx);
RigElementType eType = m_femPart->elementType(neighborElmIdx);
posKFaceIdx = RigFemTypes::oppositeFace(eType, neighborNegFaceIdx);
// Now where is posJFace of the new K cell ?
posJFaceIdx = perpendicularFaceInDirection(startElmInKNormalJ, neighborNegFaceIdx, neighborElmIdx);
posIFaceIdx = perpendicularFaceInDirection(startElmInKNormalI, neighborNegFaceIdx, neighborElmIdx);
// Step to neighbor
elmIdxInK = neighborElmIdx;
}
}
if (kCoord > m_elmentIJKCounts[2]) m_elmentIJKCounts[2] = kCoord;
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
int RigFemPartGrid::findElmIdxOfGridCornerClosestToOrigo()
{
const std::vector<int>& possibleGridCorners = m_femPart->possibleGridCornerElements();
size_t possibleCornerCount = possibleGridCorners.size();
const std::vector<cvf::Vec3f>& nodeCoordinates = m_femPart->nodes().coordinates;
int elmIdxToClosestCorner = -1;
// Find corner cell closest to origo
double minSqDistance = 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 < minSqDistance)
{
minSqDistance = distSq;
elmIdxToClosestCorner = elmIdx;
}
}
return elmIdxToClosestCorner;
}
//--------------------------------------------------------------------------------------------------
/// Find the face that is not perpFaceIdx or its opposite, and has normal closest to direction
//--------------------------------------------------------------------------------------------------
int RigFemPartGrid::perpendicularFaceInDirection(cvf::Vec3f direction, int perpFaceIdx, int elmIdx)
{
RigElementType eType = m_femPart->elementType(elmIdx);
int faceCount = RigFemTypes::elmentFaceCount(eType);
int oppFace = RigFemTypes::oppositeFace(eType, perpFaceIdx);
double minDiffSqLength = HUGE_VAL;
cvf::Vec3f faceNormal;
direction.normalize();
int bestFace = -1;
for (int faceIdx = 0; faceIdx < faceCount; ++faceIdx)
{
if (faceIdx == perpFaceIdx || faceIdx == oppFace) continue;
faceNormal = m_femPart->faceNormal(elmIdx, faceIdx);
faceNormal.normalize();
float diffSqLength = (direction - faceNormal).lengthSquared();
if (diffSqLength < minDiffSqLength)
{
bestFace = faceIdx;
minDiffSqLength = diffSqLength;
}
}
return bestFace;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigFemPartGrid::generateStructGridData()
{
// 1. Calculate neighbors for each element
// record the ones with 3 or fewer neighbors as possible grid corners
// 2. Loop over the possible corner cells,
// find the one that corresponds to IJK = 000
// by Determining what surfs correspond to NEG IJK surfaces in that element,
// and that none of those faces have a neighbor
// 4. Assign IJK = 000 to that element
// Store IJK in elm idx array
// 5. Loop along POS I surfaces increment I for each element and assign IJK
// when at end, go to POS J neighbor, increment J, repeat above.
// etc for POS Z
// Find max IJK as you go,
// also assert that there are no NEG I/NEG J/NEG Z neighbors when starting on a new row
// (Need to find min, and offset IJK values if there exists such)
// 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;
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::gridPointCountI() const
{
CVF_ASSERT(false);
return cvf::UNDEFINED_SIZE_T;
return m_elmentIJKCounts[0] + 1;
}
//--------------------------------------------------------------------------------------------------
@@ -127,8 +319,7 @@ size_t RigFemPartGrid::gridPointCountI() const
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::gridPointCountJ() const
{
CVF_ASSERT(false);
return cvf::UNDEFINED_SIZE_T;
return m_elmentIJKCounts[1] + 1;
}
//--------------------------------------------------------------------------------------------------
@@ -136,8 +327,7 @@ size_t RigFemPartGrid::gridPointCountJ() const
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::gridPointCountK() const
{
CVF_ASSERT(false);
return cvf::UNDEFINED_SIZE_T;
return m_elmentIJKCounts[2] + 1;
}
//--------------------------------------------------------------------------------------------------
@@ -190,8 +380,11 @@ size_t RigFemPartGrid::cellIndexFromIJK(size_t i, size_t j, size_t k) const
//--------------------------------------------------------------------------------------------------
bool RigFemPartGrid::ijkFromCellIndex(size_t cellIndex, size_t* i, size_t* j, size_t* k) const
{
CVF_ASSERT(false);
return false;
*i = m_ijkPrElement[cellIndex][0];
*j = m_ijkPrElement[cellIndex][1];
*k = m_ijkPrElement[cellIndex][2];
return true;
}
//--------------------------------------------------------------------------------------------------
@@ -246,5 +439,3 @@ cvf::Vec3d RigFemPartGrid::gridPointCoordinate(size_t i, size_t j, size_t k) con
CVF_ASSERT(false);
return cvf::Vec3d::ZERO;
}

12
ApplicationCode/GeoMech/GeoMechDataModel/RigFemPartGrid.h
View File

@@ -49,10 +49,16 @@ public:
private:
void generateStructGridData();
static FaceType findGridFace(cvf::Vec3d faceNormal);
void generateStructGridData();
int findElmIdxOfGridCornerClosestToOrigo();
int perpendicularFaceInDirection(cvf::Vec3f direction, int perpFaceIdx, int elmIdx);
RigFemPart* m_femPart;
std::vector<cvf::Vec3i> m_ijkPrElement;
cvf::Vec3st m_elmentIJKCounts;
RigFemPart* m_femPart;
};

2
ApplicationCode/GeoMech/GeoMechDataModel/RigFemTypes.h
View File

@@ -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);
static int oppositeFace(RigElementType elmType, int faceIdx);
};