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

Neighbor calculation and invisible face removal

Browse Source 
As a step on the road to range filtering of Geomech cases, the removal
of internal faces is implemented, and fairly optimized.
We might consider paralellization, but reading the elements is the real
bottleneck.
Memory usage should be looked at.
This commit is contained in:
Jacob Støren 2015-05-26 14:02:25 +02:00
parent f32fc130fc
commit 0d56ee060e
5 changed files with 103 additions and 58 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

93
ApplicationCode/GeoMech/GeoMechDataModel/RigFemFaceComparator.h
View File

@ -28,70 +28,65 @@ class RigFemFaceComparator
public:
void setMainFace(const int* elmNodes, const int * localFaceIndices, int faceNodeCount)
{
canonizeFace(elmNodes, localFaceIndices, faceNodeCount, &m_canonizedMainFaceIdxes);
}
bool isSameButOposite(const int* elmNodes, const int * localFaceIndices, int faceNodeCount)
{
canonizeOpositeFace(elmNodes, localFaceIndices, faceNodeCount, &m_canonizedOtherFaceIdxes);
return m_canonizedMainFaceIdxes == m_canonizedOtherFaceIdxes;
}
private:
void canonizeFace( const int* elmNodes,
const int * localFaceIndices,
int faceNodeCount,
std::vector<int>* canonizedFace)
{
canonizedFace->resize(faceNodeCount);
int minNodeIdx = INT_MAX;
int faceIdxToMinNodeIdx = 0;
m_canonizedMainFaceIdxes.resize(faceNodeCount);
m_minMainFaceNodeIdx = INT_MAX;
m_faceIdxToMinMainFaceNodeIdx = 0;
for(int fnIdx = 0; fnIdx < faceNodeCount; ++fnIdx)
{
int nodeIdx = elmNodes[localFaceIndices[fnIdx]];
(*canonizedFace)[fnIdx] = nodeIdx;
if (nodeIdx < minNodeIdx)
m_canonizedMainFaceIdxes[fnIdx] = nodeIdx;
if (nodeIdx < m_minMainFaceNodeIdx)
{
m_minMainFaceNodeIdx = nodeIdx;
m_faceIdxToMinMainFaceNodeIdx = fnIdx;
}
}
}
bool isSameButOposite(const int* elmNodes, const int * localFaceIndices, int faceNodeCount)
{
if (faceNodeCount != m_canonizedMainFaceIdxes.size()) return false;
// Find min node index in face
int minNodeIdx = INT_MAX;
int faceIdxToMinNodeIdx = 0;
for (int fnIdx = 0; fnIdx < faceNodeCount; ++fnIdx)
{
int nodeIdx = elmNodes[localFaceIndices[fnIdx]];
if (nodeIdx < minNodeIdx)
{
minNodeIdx = nodeIdx;
faceIdxToMinNodeIdx = fnIdx;
}
}
std::rotate(canonizedFace->begin(),
canonizedFace->begin() + faceIdxToMinNodeIdx,
canonizedFace->end());
}
void canonizeOpositeFace(const int* elmNodes,
const int * localFaceIndices,
int faceNodeCount,
std::vector<int>* canonizedFace)
{
canonizedFace->resize(faceNodeCount);
int minNodeIdx = INT_MAX;
int faceIdxToMinNodeIdx = 0;
int canFaceIdx = 0;
for(int fnIdx = faceNodeCount -1; fnIdx >= 0; --fnIdx, ++canFaceIdx)
// Compare faces
{
int nodeIdx = elmNodes[localFaceIndices[fnIdx]];
(*canonizedFace)[canFaceIdx] = nodeIdx;
if (nodeIdx < minNodeIdx)
{
minNodeIdx = nodeIdx;
faceIdxToMinNodeIdx = canFaceIdx;
}
}
if (minNodeIdx != m_minMainFaceNodeIdx ) return false;
std::rotate(canonizedFace->begin(),
canonizedFace->begin() + faceIdxToMinNodeIdx,
canonizedFace->end());
int canFaceIdx = m_faceIdxToMinMainFaceNodeIdx;
int fnIdx = faceIdxToMinNodeIdx;
int count = 0;
for (; count < faceNodeCount;
--fnIdx, ++canFaceIdx, ++count)
{
if (fnIdx < 0) fnIdx = faceNodeCount - 1;
if (canFaceIdx == faceNodeCount) canFaceIdx = 0;
if (elmNodes[localFaceIndices[fnIdx]] != m_canonizedMainFaceIdxes[canFaceIdx]) return false;
}
return true;
}
}
private:
std::vector<int> m_canonizedMainFaceIdxes;
std::vector<int> m_canonizedOtherFaceIdxes;
int m_minMainFaceNodeIdx;
int m_faceIdxToMinMainFaceNodeIdx;
};

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

@ -101,7 +101,7 @@ void RigFemPart::calculateNodeToElmRefs()
{
int elmNodeCount = RigFemTypes::elmentNodeCount(elementType(eIdx));
const int* elmNodes = connectivities(eIdx);
for (int enIdx = 0; enIdx < elmNodeCount; enIdx)
for (int enIdx = 0; enIdx < elmNodeCount; ++enIdx)
{
m_nodeToElmRefs[elmNodes[enIdx]].push_back(eIdx);
}
@ -144,7 +144,8 @@ void RigFemPart::calculateElmNeighbors()
// Calculate elm neighbors: elmIdxs matching each face of the element
RigFemFaceComparator fComp; // Outside loop to avoid memory alloc/dealloc. Rember to set as private in opm parallelization
std::vector<int> candidates;//
m_elmNeighbors.resize(this->elementCount());
for (cvf::uint eIdx = 0; eIdx < this->elementCount(); ++eIdx)
@ -162,18 +163,51 @@ void RigFemPart::calculateElmNeighbors()
const int* localFaceIndices = RigFemTypes::localElmNodeIndicesForFace(elmType, faceIdx, &faceNodeCount);
// Get neighbor candidates
int firstNodeIdxOfFace = elmNodes[localFaceIndices[0]];
int neigborCandidateCount = this->numElementsUsingNode(firstNodeIdxOfFace);
const size_t* candidates = this->elementsUsingNode(firstNodeIdxOfFace);
if (neigborCandidateCount)
candidates.clear();
{
int firstNodeIdxOfFace = elmNodes[localFaceIndices[0]];
int candidateCount1 = this->numElementsUsingNode(firstNodeIdxOfFace);
const size_t* candidates1 = this->elementsUsingNode(firstNodeIdxOfFace);
if (candidateCount1)
{
// Get neighbor candidates from the diagonal node
int thirdNodeIdxOfFace = elmNodes[localFaceIndices[3]];
int candidateCount2 = this->numElementsUsingNode(thirdNodeIdxOfFace);
const size_t* 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
int idx1 = 0;
int idx2 = 0;
while (idx1 < candidateCount1 && idx2 < candidateCount2)
{
if (candidates1[idx1] < candidates2[idx2]){ ++idx1; continue; }
if (candidates1[idx1] > candidates2[idx2]){ ++idx2; continue; }
if (candidates1[idx1] == candidates2[idx2])
{
if (candidates1[idx1] != eIdx)
{
candidates.push_back(candidates1[idx1]);
}
++idx1; ++idx2;
}
}
}
}
if (candidates.size())
{
fComp.setMainFace(elmNodes, localFaceIndices, faceNodeCount);
}
// Check if any of the neighbor candidates faces matches
for (int nbcIdx = 0; nbcIdx < neigborCandidateCount; ++nbcIdx)
for (int nbcIdx = 0; nbcIdx < candidates.size(); ++nbcIdx)
{
int nbcElmIdx = candidates[neigborCandidateCount];
int nbcElmIdx = candidates[nbcIdx];
RigElementType nbcElmType = this->elementType(nbcElmIdx);
const int* nbcElmNodes = this->connectivities(nbcElmIdx);

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

@ -66,7 +66,7 @@ public:
int numElementsUsingNode(int nodeIndex);
void assertElmNeighborsIsCalculated();
int elementNeighbor(int elementIndex, int faceIndex)
int elementNeighbor(int elementIndex, int faceIndex) const
{ return m_elmNeighbors[elementIndex].idxToNeighborElmPrFace[faceIndex]; }
const RigFemPartGrid* structGrid();

7
ApplicationCode/GeoMech/GeoMechDataModel/RigGeoMechCaseData.cpp
View File

@ -87,6 +87,13 @@ bool RigGeoMechCaseData::openAndReadFemParts()
m_femPartResults[pIdx] = new RigFemPartResults;
m_femPartResults[pIdx]->initResultStages(stepNames);
}
// Calculate derived Fem data
for (int pIdx = 0; pIdx < m_femParts->partCount(); ++pIdx)
{
m_femParts->part(pIdx)->assertNodeToElmIndicesIsCalculated();
m_femParts->part(pIdx)->assertElmNeighborsIsCalculated();
}
return true;
}
}

9
ApplicationCode/GeoMech/GeoMechVisualization/RivFemPartGeometryGenerator.cpp
View File

@ -162,6 +162,8 @@ void RivFemPartGeometryGenerator::computeArrays()
m_quadVerticesToNodeIdx.reserve(estimatedQuadVxCount);
m_quadVerticesToGlobalElmNodeIdx.reserve(estimatedQuadVxCount);
cvf::Vec3d offset = Vec3d::ZERO; //m_part->displayModelOffset();
const std::vector<cvf::Vec3f>& nodeCoordinates = m_part->nodes().coordinates;
@ -178,6 +180,13 @@ void RivFemPartGeometryGenerator::computeArrays()
for (int lfIdx = 0; lfIdx < faceCount; ++lfIdx)
{
int elmNeighbor = m_part->elementNeighbor(elmIdx, lfIdx);
if (elmNeighbor != -1 && (m_elmVisibility.isNull() || (*m_elmVisibility)[elmNeighbor]))
{
continue; // Invisible face
}
int faceNodeCount = 0;
const int* localElmNodeIndicesForFace = RigFemTypes::localElmNodeIndicesForFace(eType, lfIdx, &faceNodeCount);
if (faceNodeCount == 4)