Prep Geom range: Neighbor calculation

A checkpoint before starting to use the neighbor calculations
This commit is contained in:
Jacob Støren 2015-05-26 08:57:53 +02:00
parent b280196b7c
commit f32fc130fc
7 changed files with 468 additions and 157 deletions

View File

@ -24,6 +24,9 @@ add_library( ${PROJECT_NAME}
RigFemScalarResultFrames.cpp
RigFemNativeStatCalc.h
RigFemNativeStatCalc.cpp
RigFemFaceComparator.h
RigFemPartGrid.h
RigFemPartGrid.cpp
)
target_link_libraries( ${PROJECT_NAME} LibCore ResultStatisticsCache)

View File

@ -0,0 +1,97 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2015- Statoil ASA
// Copyright (C) 2015- Ceetron Solutions AS
//
// ResInsight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE.
//
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#pragma once
#include <vector>
#include <algorithm>
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;
for(int fnIdx = 0; fnIdx < faceNodeCount; ++fnIdx)
{
int nodeIdx = elmNodes[localFaceIndices[fnIdx]];
(*canonizedFace)[fnIdx] = nodeIdx;
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)
{
int nodeIdx = elmNodes[localFaceIndices[fnIdx]];
(*canonizedFace)[canFaceIdx] = nodeIdx;
if (nodeIdx < minNodeIdx)
{
minNodeIdx = nodeIdx;
faceIdxToMinNodeIdx = canFaceIdx;
}
}
std::rotate(canonizedFace->begin(),
canonizedFace->begin() + faceIdxToMinNodeIdx,
canonizedFace->end());
}
std::vector<int> m_canonizedMainFaceIdxes;
std::vector<int> m_canonizedOtherFaceIdxes;
};

View File

@ -19,6 +19,7 @@
#include "RigFemPart.h"
#include "RigFemPartGrid.h"
//--------------------------------------------------------------------------------------------------
///
@ -78,162 +79,124 @@ const RigFemPartGrid* RigFemPart::structGrid()
return m_structGrid.p();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigFemPartGrid::RigFemPartGrid(RigFemPart* femPart)
void RigFemPart::assertNodeToElmIndicesIsCalculated()
{
m_femPart = femPart;
generateStructGridData();
if (m_nodeToElmRefs.size() != nodes().nodeIds.size())
{
this->calculateNodeToElmRefs();
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigFemPartGrid::~RigFemPartGrid()
void RigFemPart::calculateNodeToElmRefs()
{
m_nodeToElmRefs.resize(nodes().nodeIds.size());
for (size_t eIdx = 0; eIdx < m_elementId.size(); ++eIdx)
{
int elmNodeCount = RigFemTypes::elmentNodeCount(elementType(eIdx));
const int* elmNodes = connectivities(eIdx);
for (int enIdx = 0; enIdx < elmNodeCount; enIdx)
{
m_nodeToElmRefs[elmNodes[enIdx]].push_back(eIdx);
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigFemPartGrid::generateStructGridData()
const size_t* RigFemPart::elementsUsingNode(int nodeIndex)
{
return &(m_nodeToElmRefs[nodeIndex][0]);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::gridPointCountI() const
int RigFemPart::numElementsUsingNode(int nodeIndex)
{
CVF_ASSERT(false);
return cvf::UNDEFINED_SIZE_T;
return static_cast<int>(m_nodeToElmRefs[nodeIndex].size());
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::gridPointCountJ() const
void RigFemPart::assertElmNeighborsIsCalculated()
{
CVF_ASSERT(false);
return cvf::UNDEFINED_SIZE_T;
if (m_elmNeighbors.size() != m_elementId.size())
{
this->calculateElmNeighbors();
}
}
#include "RigFemFaceComparator.h"
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::gridPointCountK() const
void RigFemPart::calculateElmNeighbors()
{
CVF_ASSERT(false);
return cvf::UNDEFINED_SIZE_T;
// 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
m_elmNeighbors.resize(this->elementCount());
for (cvf::uint eIdx = 0; eIdx < this->elementCount(); ++eIdx)
{
RigElementType elmType = this->elementType(eIdx);
const int* elmNodes = this->connectivities(eIdx);
int faceCount = RigFemTypes::elmentFaceCount(elmType);
for (int faceIdx = 0; faceIdx < faceCount; ++faceIdx)
{
m_elmNeighbors[eIdx].idxToNeighborElmPrFace[faceIdx] = -1;
int faceNodeCount = 0;
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)
{
fComp.setMainFace(elmNodes, localFaceIndices, faceNodeCount);
}
// Check if any of the neighbor candidates faces matches
for (int nbcIdx = 0; nbcIdx < neigborCandidateCount; ++nbcIdx)
{
int nbcElmIdx = candidates[neigborCandidateCount];
RigElementType nbcElmType = this->elementType(nbcElmIdx);
const int* nbcElmNodes = this->connectivities(nbcElmIdx);
int nbcFaceCount = RigFemTypes::elmentFaceCount(nbcElmType);
bool isNeighborFound = false;
for (int nbcFaceIdx = 0; nbcFaceIdx < nbcFaceCount; ++nbcFaceIdx)
{
int nbcFaceNodeCount = 0;
const int* nbcLocalFaceIndices = RigFemTypes::localElmNodeIndicesForFace(nbcElmType, nbcFaceIdx, &nbcFaceNodeCount);
// Compare faces
if (fComp.isSameButOposite(nbcElmNodes, nbcLocalFaceIndices, nbcFaceNodeCount))
{
m_elmNeighbors[eIdx].idxToNeighborElmPrFace[faceIdx] = nbcElmIdx;
isNeighborFound = true;
break;
}
}
if (isNeighborFound) break;
}
}
}
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartGrid::isCellValid(size_t i, size_t j, size_t k) const
{
CVF_ASSERT(false);
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigFemPartGrid::minCoordinate() const
{
CVF_ASSERT(false);
return cvf::Vec3d::ZERO;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigFemPartGrid::maxCoordinate() const
{
CVF_ASSERT(false);
return cvf::Vec3d::ZERO;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartGrid::cellIJKNeighbor(size_t i, size_t j, size_t k, FaceType face, size_t* neighborCellIndex) const
{
CVF_ASSERT(false);
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::cellIndexFromIJK(size_t i, size_t j, size_t k) const
{
CVF_ASSERT(false);
return cvf::UNDEFINED_SIZE_T;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartGrid::ijkFromCellIndex(size_t cellIndex, size_t* i, size_t* j, size_t* k) const
{
CVF_ASSERT(false);
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartGrid::cellIJKFromCoordinate(const cvf::Vec3d& coord, size_t* i, size_t* j, size_t* k) const
{
CVF_ASSERT(false);
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigFemPartGrid::cellCornerVertices(size_t cellIndex, cvf::Vec3d vertices[8]) const
{
CVF_ASSERT(false);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigFemPartGrid::cellCentroid(size_t cellIndex) const
{
CVF_ASSERT(false);
return cvf::Vec3d::ZERO;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigFemPartGrid::cellMinMaxCordinates(size_t cellIndex, cvf::Vec3d* minCoordinate, cvf::Vec3d* maxCoordinate) const
{
CVF_ASSERT(false);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::gridPointIndexFromIJK(size_t i, size_t j, size_t k) const
{
CVF_ASSERT(false);
return cvf::UNDEFINED_SIZE_T;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigFemPartGrid::gridPointCoordinate(size_t i, size_t j, size_t k) const
{
CVF_ASSERT(false);
return cvf::Vec3d::ZERO;
}

View File

@ -35,6 +35,8 @@ class RigFemPartNodes
public:
std::vector<int> nodeIds;
std::vector<cvf::Vec3f> coordinates;
};
class RigFemPart : public cvf::Object
@ -59,10 +61,18 @@ public:
RigFemPartNodes& nodes() {return m_nodes;}
const RigFemPartNodes& nodes() const {return m_nodes;}
const RigFemPartGrid* structGrid();
void assertNodeToElmIndicesIsCalculated();
const size_t* elementsUsingNode(int nodeIndex);
int numElementsUsingNode(int nodeIndex);
void assertElmNeighborsIsCalculated();
int elementNeighbor(int elementIndex, int faceIndex)
{ return m_elmNeighbors[elementIndex].idxToNeighborElmPrFace[faceIndex]; }
const RigFemPartGrid* structGrid();
private:
int m_elementPartId;
std::vector<int> m_elementId;
std::vector<RigElementType> m_elementTypes;
std::vector<size_t> m_elementConnectivityStartIndices;
@ -71,37 +81,11 @@ private:
RigFemPartNodes m_nodes;
cvf::ref<RigFemPartGrid> m_structGrid;
void calculateNodeToElmRefs();
std::vector<std::vector<size_t>> m_nodeToElmRefs; // Needs a more memory friendly structure
void calculateElmNeighbors();
struct Neighbors { int idxToNeighborElmPrFace[6]; };
std::vector< Neighbors > m_elmNeighbors;
};
#include "cvfStructGrid.h"
class RigFemPartGrid : public cvf::StructGridInterface
{
public:
RigFemPartGrid(RigFemPart* femPart);
virtual ~RigFemPartGrid();
virtual size_t gridPointCountI() const;
virtual size_t gridPointCountJ() const;
virtual size_t gridPointCountK() const;
virtual bool isCellValid(size_t i, size_t j, size_t k) const;
virtual cvf::Vec3d minCoordinate() const;
virtual cvf::Vec3d maxCoordinate() const;
virtual bool cellIJKNeighbor(size_t i, size_t j, size_t k, FaceType face, size_t* neighborCellIndex) const;
virtual size_t cellIndexFromIJK(size_t i, size_t j, size_t k) const;
virtual bool ijkFromCellIndex(size_t cellIndex, size_t* i, size_t* j, size_t* k) const;
virtual bool cellIJKFromCoordinate(const cvf::Vec3d& coord, size_t* i, size_t* j, size_t* k) const;
virtual void cellCornerVertices(size_t cellIndex, cvf::Vec3d vertices[8]) const;
virtual cvf::Vec3d cellCentroid(size_t cellIndex) const;
virtual void cellMinMaxCordinates(size_t cellIndex, cvf::Vec3d* minCoordinate, cvf::Vec3d* maxCoordinate) const;
virtual size_t gridPointIndexFromIJK(size_t i, size_t j, size_t k) const;
virtual cvf::Vec3d gridPointCoordinate(size_t i, size_t j, size_t k) const;
private:
void generateStructGridData();
RigFemPart* m_femPart;
};

View File

@ -0,0 +1,203 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2015- Statoil ASA
// Copyright (C) 2015- Ceetron Solutions AS
//
// ResInsight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE.
//
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#include "RigFemPartGrid.h"
#include "RigFemPart.h"
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigFemPartGrid::RigFemPartGrid(RigFemPart* femPart)
{
m_femPart = femPart;
generateStructGridData();
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
RigFemPartGrid::~RigFemPartGrid()
{
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
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
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::gridPointCountI() const
{
CVF_ASSERT(false);
return cvf::UNDEFINED_SIZE_T;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::gridPointCountJ() const
{
CVF_ASSERT(false);
return cvf::UNDEFINED_SIZE_T;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::gridPointCountK() const
{
CVF_ASSERT(false);
return cvf::UNDEFINED_SIZE_T;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartGrid::isCellValid(size_t i, size_t j, size_t k) const
{
CVF_ASSERT(false);
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigFemPartGrid::minCoordinate() const
{
CVF_ASSERT(false);
return cvf::Vec3d::ZERO;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigFemPartGrid::maxCoordinate() const
{
CVF_ASSERT(false);
return cvf::Vec3d::ZERO;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartGrid::cellIJKNeighbor(size_t i, size_t j, size_t k, FaceType face, size_t* neighborCellIndex) const
{
CVF_ASSERT(false);
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::cellIndexFromIJK(size_t i, size_t j, size_t k) const
{
CVF_ASSERT(false);
return cvf::UNDEFINED_SIZE_T;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartGrid::ijkFromCellIndex(size_t cellIndex, size_t* i, size_t* j, size_t* k) const
{
CVF_ASSERT(false);
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
bool RigFemPartGrid::cellIJKFromCoordinate(const cvf::Vec3d& coord, size_t* i, size_t* j, size_t* k) const
{
CVF_ASSERT(false);
return false;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigFemPartGrid::cellCornerVertices(size_t cellIndex, cvf::Vec3d vertices[8]) const
{
CVF_ASSERT(false);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigFemPartGrid::cellCentroid(size_t cellIndex) const
{
CVF_ASSERT(false);
return cvf::Vec3d::ZERO;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
void RigFemPartGrid::cellMinMaxCordinates(size_t cellIndex, cvf::Vec3d* minCoordinate, cvf::Vec3d* maxCoordinate) const
{
CVF_ASSERT(false);
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
size_t RigFemPartGrid::gridPointIndexFromIJK(size_t i, size_t j, size_t k) const
{
CVF_ASSERT(false);
return cvf::UNDEFINED_SIZE_T;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::Vec3d RigFemPartGrid::gridPointCoordinate(size_t i, size_t j, size_t k) const
{
CVF_ASSERT(false);
return cvf::Vec3d::ZERO;
}

View File

@ -0,0 +1,60 @@
/////////////////////////////////////////////////////////////////////////////////
//
// Copyright (C) 2015- Statoil ASA
// Copyright (C) 2015- Ceetron Solutions AS
//
// ResInsight is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// ResInsight is distributed in the hope that it will be useful, but WITHOUT ANY
// WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE.
//
// See the GNU General Public License at <http://www.gnu.org/licenses/gpl.html>
// for more details.
//
/////////////////////////////////////////////////////////////////////////////////
#pragma once
#include "cvfStructGrid.h"
class RigFemPart;
class RigFemPartGrid : public cvf::StructGridInterface
{
public:
RigFemPartGrid(RigFemPart* femPart);
virtual ~RigFemPartGrid();
virtual size_t gridPointCountI() const;
virtual size_t gridPointCountJ() const;
virtual size_t gridPointCountK() const;
virtual bool isCellValid(size_t i, size_t j, size_t k) const;
virtual cvf::Vec3d minCoordinate() const;
virtual cvf::Vec3d maxCoordinate() const;
virtual bool cellIJKNeighbor(size_t i, size_t j, size_t k, FaceType face, size_t* neighborCellIndex) const;
virtual size_t cellIndexFromIJK(size_t i, size_t j, size_t k) const;
virtual bool ijkFromCellIndex(size_t cellIndex, size_t* i, size_t* j, size_t* k) const;
virtual bool cellIJKFromCoordinate(const cvf::Vec3d& coord, size_t* i, size_t* j, size_t* k) const;
virtual void cellCornerVertices(size_t cellIndex, cvf::Vec3d vertices[8]) const;
virtual cvf::Vec3d cellCentroid(size_t cellIndex) const;
virtual void cellMinMaxCordinates(size_t cellIndex, cvf::Vec3d* minCoordinate, cvf::Vec3d* maxCoordinate) const;
virtual size_t gridPointIndexFromIJK(size_t i, size_t j, size_t k) const;
virtual cvf::Vec3d gridPointCoordinate(size_t i, size_t j, size_t k) const;
private:
void generateStructGridData();
RigFemPart* m_femPart;
};

View File

@ -28,7 +28,8 @@
#include "RimGeoMechView.h"
#include "RimGeoMechCase.h"
#include "RigGeomechCaseData.h"
#include "RigFemPart.h"
#include "RigFemPartGrid.h"
CAF_PDM_SOURCE_INIT(RimCellRangeFilterCollection, "CellRangeFilterCollection");