mirror of
https://github.com/OPM/ResInsight.git
synced 2025-02-25 18:55:39 -06:00
#815 First working interpolation scheme, so values are not insane
This commit is contained in:
@@ -92,6 +92,64 @@ class RivIntersectionVertexWeights
|
||||
public:
|
||||
explicit RivIntersectionVertexWeights(): m_count(0) {}
|
||||
|
||||
/*
|
||||
v111 k11 v11 v112 v211 k21 v21 v212
|
||||
+------+-----+ +--------+------+
|
||||
| |
|
||||
|k1 |k2
|
||||
| k |
|
||||
v1+--------+------------+v2
|
||||
| |
|
||||
| |
|
||||
| |
|
||||
+------+----+ +--------+-------+
|
||||
v121 k12 v12 v122 v221 k22 v22 v222
|
||||
|
||||
Where the k's are normalized distances along the edge, from the edge start vertex .
|
||||
|
||||
This is the interpolation sceme:
|
||||
|
||||
v = (1 - k )* v1 + k *v2;
|
||||
|
||||
v1 = (1 - k1 )* v11 + k1*v12;
|
||||
v2 = (1 - k2 )* v21 + k2*v22;
|
||||
|
||||
v11 = (1 - k11)* v111 + k11*v112;
|
||||
v12 = (1 - k12)* v121 + k12*v122;
|
||||
v21 = (1 - k21)* v211 + k21*v212;
|
||||
v22 = (1 - k22)* v221 + k22*v222;
|
||||
|
||||
Substitution and reorganizing gives the expressions seen below.
|
||||
|
||||
*/
|
||||
|
||||
explicit RivIntersectionVertexWeights( size_t edgeVx111, size_t edgeVx112, double k11,
|
||||
size_t edgeVx121, size_t edgeVx122, double k12,
|
||||
size_t edgeVx211, size_t edgeVx212, double k21,
|
||||
size_t edgeVx221, size_t edgeVx222, double k22,
|
||||
double k1,
|
||||
double k2,
|
||||
double k) : m_count(8)
|
||||
{
|
||||
m_vxIds[0] = (edgeVx111);
|
||||
m_vxIds[1] = (edgeVx112);
|
||||
m_vxIds[2] = (edgeVx121);
|
||||
m_vxIds[3] = (edgeVx122);
|
||||
m_vxIds[4] = (edgeVx211);
|
||||
m_vxIds[5] = (edgeVx212);
|
||||
m_vxIds[6] = (edgeVx221);
|
||||
m_vxIds[7] = (edgeVx222);
|
||||
|
||||
m_weights[7] = ((float)(k * k2 * k22));
|
||||
m_weights[6] = ((float)(k * k2 - k * k2 * k22));
|
||||
m_weights[5] = ((float)(( k - k * k2 ) * k21));
|
||||
m_weights[4] = ((float)((k * k2 - k ) * k21 - k * k2 + k));
|
||||
m_weights[3] = ((float)((1-k) * k1 * k12));
|
||||
m_weights[2] = ((float)((k-1) * k1 * k12 + ( 1 - k ) * k1));
|
||||
m_weights[1] = ((float)(( (k-1) * k1 - k + 1 ) * k11));
|
||||
m_weights[0] = ((float)(( (1-k) * k1 + k - 1 ) * k11 + ( k - 1 ) * k1 - k + 1));
|
||||
}
|
||||
|
||||
explicit RivIntersectionVertexWeights(size_t edge1Vx1, size_t edge1Vx2, double normDistFromE1V1,
|
||||
size_t edge2Vx1, size_t edge2Vx2, double normDistFromE2V1,
|
||||
double normDistFromE1Cut) : m_count(4)
|
||||
@@ -122,7 +180,7 @@ public:
|
||||
|
||||
private:
|
||||
|
||||
size_t m_vxIds[4];
|
||||
float m_weights[4];
|
||||
size_t m_vxIds[8];
|
||||
float m_weights[8];
|
||||
int m_count;
|
||||
};
|
||||
|
||||
@@ -256,7 +256,7 @@ void RivIntersectionBoxGeometryGenerator::calculateArrays()
|
||||
|
||||
for (cvf::Vec3d& corner : faceCorners) sectionBBox.add(corner);
|
||||
|
||||
// Nearly same code as IntersectionGenerator :
|
||||
// Similar code as IntersectionGenerator :
|
||||
|
||||
std::vector<size_t> columnCellCandidates;
|
||||
m_hexGrid->findIntersectingCells(sectionBBox, &columnCellCandidates);
|
||||
@@ -289,12 +289,14 @@ void RivIntersectionBoxGeometryGenerator::calculateArrays()
|
||||
caf::HexGridIntersectionTools::clipTrianglesBetweenTwoParallelPlanes(hexPlaneCutTriangleVxes, isTriangleEdgeCellContour, p1Plane, p2Plane,
|
||||
&clippedTriangleVxes_once, &isClippedTriEdgeCellContour_once);
|
||||
|
||||
for (caf::HexGridIntersectionTools::ClipVx& clvx : clippedTriangleVxes_once) if (!clvx.isVxIdsNative) clvx.derivedVxLevel = 0;
|
||||
|
||||
std::vector<caf::HexGridIntersectionTools::ClipVx> clippedTriangleVxes;
|
||||
std::vector<bool> isClippedTriEdgeCellContour;
|
||||
|
||||
caf::HexGridIntersectionTools::clipTrianglesBetweenTwoParallelPlanes(clippedTriangleVxes_once, isClippedTriEdgeCellContour_once, p3Plane, p4Plane,
|
||||
&clippedTriangleVxes, &isClippedTriEdgeCellContour);
|
||||
|
||||
for (caf::HexGridIntersectionTools::ClipVx& clvx : clippedTriangleVxes) if (!clvx.isVxIdsNative && clvx.derivedVxLevel == -1) clvx.derivedVxLevel = 1;
|
||||
|
||||
size_t clippedTriangleCount = clippedTriangleVxes.size()/3;
|
||||
|
||||
@@ -346,19 +348,92 @@ void RivIntersectionBoxGeometryGenerator::calculateArrays()
|
||||
}
|
||||
else
|
||||
{
|
||||
#if 0
|
||||
caf::HexGridIntersectionTools::ClipVx cvx1 = hexPlaneCutTriangleVxes[cvx.clippedEdgeVx1Id];
|
||||
caf::HexGridIntersectionTools::ClipVx cvx2 = hexPlaneCutTriangleVxes[cvx.clippedEdgeVx2Id];
|
||||
caf::HexGridIntersectionTools::ClipVx cvx1;
|
||||
caf::HexGridIntersectionTools::ClipVx cvx2;
|
||||
|
||||
m_triVxToCellCornerWeights.push_back(
|
||||
RivIntersectionVertexWeights(cvx1.clippedEdgeVx1Id, cvx1.clippedEdgeVx2Id, cvx1.normDistFromEdgeVx1,
|
||||
cvx2.clippedEdgeVx1Id, cvx2.clippedEdgeVx2Id, cvx2.normDistFromEdgeVx1,
|
||||
cvx.normDistFromEdgeVx1));
|
||||
#else
|
||||
m_triVxToCellCornerWeights.push_back(RivIntersectionVertexWeights());
|
||||
#endif
|
||||
if (cvx.derivedVxLevel == 0)
|
||||
{
|
||||
cvx1 = hexPlaneCutTriangleVxes[cvx.clippedEdgeVx1Id];
|
||||
cvx2 = hexPlaneCutTriangleVxes[cvx.clippedEdgeVx2Id];
|
||||
}
|
||||
else if(cvx.derivedVxLevel == 1)
|
||||
{
|
||||
cvx1 = clippedTriangleVxes_once[cvx.clippedEdgeVx1Id];
|
||||
cvx2 = clippedTriangleVxes_once[cvx.clippedEdgeVx2Id];
|
||||
}
|
||||
else
|
||||
{
|
||||
CVF_ASSERT(false);
|
||||
}
|
||||
|
||||
if(cvx1.isVxIdsNative && cvx2.isVxIdsNative)
|
||||
{
|
||||
m_triVxToCellCornerWeights.push_back(
|
||||
RivIntersectionVertexWeights(cvx1.clippedEdgeVx1Id, cvx1.clippedEdgeVx2Id, cvx1.normDistFromEdgeVx1,
|
||||
cvx2.clippedEdgeVx1Id, cvx2.clippedEdgeVx2Id, cvx2.normDistFromEdgeVx1,
|
||||
cvx.normDistFromEdgeVx1));
|
||||
}
|
||||
else
|
||||
{
|
||||
caf::HexGridIntersectionTools::ClipVx cvx11;
|
||||
caf::HexGridIntersectionTools::ClipVx cvx12;
|
||||
caf::HexGridIntersectionTools::ClipVx cvx21;
|
||||
caf::HexGridIntersectionTools::ClipVx cvx22;
|
||||
|
||||
if(cvx1.isVxIdsNative)
|
||||
{
|
||||
cvx11 = cvx1;
|
||||
cvx12 = cvx1;
|
||||
}
|
||||
else if(cvx1.derivedVxLevel == 0)
|
||||
{
|
||||
cvx11 = hexPlaneCutTriangleVxes[cvx1.clippedEdgeVx1Id];
|
||||
cvx12 = hexPlaneCutTriangleVxes[cvx1.clippedEdgeVx2Id];
|
||||
}
|
||||
else if(cvx2.derivedVxLevel == 1)
|
||||
{
|
||||
cvx11 = clippedTriangleVxes_once[cvx1.clippedEdgeVx1Id];
|
||||
cvx12 = clippedTriangleVxes_once[cvx1.clippedEdgeVx2Id];
|
||||
}
|
||||
else
|
||||
{
|
||||
CVF_ASSERT(false);
|
||||
}
|
||||
|
||||
|
||||
if(cvx2.isVxIdsNative)
|
||||
{
|
||||
cvx21 = cvx2;
|
||||
cvx22 = cvx2;
|
||||
}
|
||||
else if(cvx2.derivedVxLevel == 0)
|
||||
{
|
||||
cvx21 = hexPlaneCutTriangleVxes[cvx2.clippedEdgeVx1Id];
|
||||
cvx22 = hexPlaneCutTriangleVxes[cvx2.clippedEdgeVx2Id];
|
||||
}
|
||||
else if(cvx2.derivedVxLevel == 1)
|
||||
{
|
||||
cvx21 = clippedTriangleVxes_once[cvx2.clippedEdgeVx1Id];
|
||||
cvx22 = clippedTriangleVxes_once[cvx2.clippedEdgeVx2Id];
|
||||
|
||||
}
|
||||
else
|
||||
{
|
||||
CVF_ASSERT(false);
|
||||
}
|
||||
|
||||
CVF_TIGHT_ASSERT(cvx11.isVxIdsNative && cvx12.isVxIdsNative && cvx21.isVxIdsNative && cvx22.isVxIdsNative);
|
||||
|
||||
m_triVxToCellCornerWeights.push_back(
|
||||
RivIntersectionVertexWeights(cvx11.clippedEdgeVx1Id, cvx11.clippedEdgeVx2Id, cvx11.normDistFromEdgeVx1,
|
||||
cvx12.clippedEdgeVx1Id, cvx12.clippedEdgeVx2Id, cvx2.normDistFromEdgeVx1,
|
||||
cvx21.clippedEdgeVx1Id, cvx21.clippedEdgeVx2Id, cvx21.normDistFromEdgeVx1,
|
||||
cvx22.clippedEdgeVx1Id, cvx22.clippedEdgeVx2Id, cvx22.normDistFromEdgeVx1,
|
||||
cvx1.normDistFromEdgeVx1,
|
||||
cvx2.normDistFromEdgeVx1,
|
||||
cvx.normDistFromEdgeVx1));
|
||||
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user