#5437 Use marching tets instead of marching cubes,

to make the cell border edges consistent
This commit is contained in:
Jacob Støren 2020-01-29 15:44:36 +01:00
parent f5b59a6bd9
commit ddb5d0fd30
3 changed files with 392 additions and 24 deletions

View File

@ -190,11 +190,12 @@ void RivSurfaceIntersectionGeometryGenerator::calculateArrays()
m_hexGrid->cellCornerIndices( globalCellIdx, &cornerIndices[0] );
hexPlaneCutTriangleVxes.clear();
int triangleCount = caf::HexGridIntersectionTools::planeHexIntersectionMC( plane,
&cellCorners[0],
&cornerIndices[0],
&hexPlaneCutTriangleVxes,
&cellFaceForEachTriangleEdge );
cellFaceForEachTriangleEdge.clear();
int triangleCount = caf::HexGridIntersectionTools::planeHexIntersectionMCTet( plane,
&cellCorners[0],
&cornerIndices[0],
&hexPlaneCutTriangleVxes,
&cellFaceForEachTriangleEdge );
if ( triangleCount == 0 ) continue;

View File

@ -540,26 +540,26 @@ void HexGridIntersectionTools::clipPlanarTrianglesWithInPlaneTriangle(const std:
// Creating a plane for each of the edges of the clipping triangle
std::array<cvf::Plane, 3> clipTrianglePlanes;
clipTrianglePlanes[0] = createPlaneFromEdgeAndPointInNormalDirection ( tp1, tp2, tp3 );
clipTrianglePlanes[1] = createPlaneFromEdgeAndPointInNormalDirection ( tp2, tp3, tp1 );
clipTrianglePlanes[2] = createPlaneFromEdgeAndPointInNormalDirection ( tp3, tp1, tp2 );
clipTrianglePlanes[0] = createPlaneFromEdgeAndPointInNormalDirection ( tp1, tp2, tp3 );
clipTrianglePlanes[1] = createPlaneFromEdgeAndPointInNormalDirection ( tp2, tp3, tp1 );
clipTrianglePlanes[2] = createPlaneFromEdgeAndPointInNormalDirection ( tp3, tp1, tp2 );
#define reserveSize 60
#define reserveSize 60
std::vector<cvf::Vec3d> currentInputTriangleVxes;
currentInputTriangleVxes.reserve(reserveSize);
std::vector<int> currentInputCellFaceForEachTriangleEdge;
currentInputCellFaceForEachTriangleEdge.reserve(reserveSize);
std::vector<cvf::Vec3d> currentOutputTriangleVxes;
currentOutputTriangleVxes.reserve(reserveSize);
std::vector<int> currentOutputCellFaceForEachTriangleEdge;
currentOutputCellFaceForEachTriangleEdge.reserve(reserveSize);
std::vector<cvf::Vec3d> currentInputTriangleVxes;
currentInputTriangleVxes.reserve(reserveSize);
std::vector<int> currentInputCellFaceForEachTriangleEdge;
currentInputCellFaceForEachTriangleEdge.reserve(reserveSize);
std::vector<cvf::Vec3d> currentOutputTriangleVxes;
currentOutputTriangleVxes.reserve(reserveSize);
std::vector<int> currentOutputCellFaceForEachTriangleEdge;
currentOutputCellFaceForEachTriangleEdge.reserve(reserveSize);
for( size_t tIdx = 0; tIdx < triangleCount; ++tIdx )
for( size_t tIdx = 0; tIdx < triangleCount; ++tIdx )
{
size_t triVxIdx = tIdx * 3;
currentInputTriangleVxes.clear();
currentInputTriangleVxes.clear();
currentInputCellFaceForEachTriangleEdge.clear();
currentOutputTriangleVxes.clear();
currentOutputCellFaceForEachTriangleEdge.clear();
@ -572,14 +572,14 @@ void HexGridIntersectionTools::clipPlanarTrianglesWithInPlaneTriangle(const std:
currentOutputCellFaceForEachTriangleEdge.push_back(cellFaceForEachTriangleEdge[triVxIdx + 1]);
currentOutputCellFaceForEachTriangleEdge.push_back(cellFaceForEachTriangleEdge[triVxIdx + 2]);
ClipVx newVx1;
newVx1.isVxIdsNative = false;
ClipVx newVx2;
newVx2.isVxIdsNative = false;
ClipVx newVx1;
newVx1.isVxIdsNative = false;
ClipVx newVx2;
newVx2.isVxIdsNative = false;
for ( int planeIdx = 0; planeIdx < 3; ++planeIdx )
{
currentInputTriangleVxes.swap(currentOutputTriangleVxes);
currentInputTriangleVxes.swap(currentOutputTriangleVxes);
currentInputCellFaceForEachTriangleEdge.swap(currentOutputCellFaceForEachTriangleEdge);
currentOutputTriangleVxes.clear();
@ -1203,5 +1203,360 @@ int HexGridIntersectionTools::planeHexIntersectionMC(const cvf::Plane& plane,
//--------------------------------------------------------------------------------------------------
/// Based on description and implementation from Paul Bourke:
///
/// http://paulbourke.net/geometry/polygonise/
///
/// Note that the element is turned inside-out compared to what we use elsewhere in caf/ResInsight
/// So the winding of all the sides are opposite.
/// 4-----4------5
/// /| /| k POS_I = 0
/// 7 8 5 9 | NEG_I = 1
/// / | / | | POS_J = 2
/// 7------6-----6 | | NEG_J = 3
/// | 0-----0--|---1 *------i POS_K = 4
/// 11 / 10 / / NEG_K = 5
/// | 3 | 1 / NO_FACE = 6
/// |/ |/ j
/// 3------2-----2
///
// The cellFaceForEachTriangleEdge refer to the edge after the corresponding triangle vertex.
/*
Based on description and implementation from Paul Bourke:
http://local.wasp.uwa.edu.au/~pbourke/geometry/polygonise/
Polygonise a tetrahedron given its vertices within a cube
This is an alternative algorithm to polygonisegrid.
It results in a smoother surface but more triangular facets.
+ 0 + 0
/|\ /|\
/ | \ / | \
/ | \ / | \
/ | \ / | \
/ | \ / 2 | \
/ | \ / __--+_ \
+-------------+ 1 3 +--__ --_ \
3 \ | / ---__ --_\
\ | / ---__-\
\ | / --+ 1
\ | /
\ | /
\|/ 2 is behind 1 and 3
+ 2
Build six tets from a cube to make sure the split direction is equal for opposite sides.
Surface normals are pointing outward.
See following comment is taken from http://www.iue.tuwien.ac.at/phd/wessner/node32.html
The decompositions of a cube into five tetrahedra yields an orientation switch of two opposite diagonal face edges of the cube.
Due to this fact, the tessellation of one cube, as part of a larger cubic grid, forces a particular tessellation of all neighboring cubes to guarantee a conformal mesh.
This means that, if such five-decompositions cubes are stacked together to a chain, the mesh of each cube must be rotated by an angle of 90 deg
The tessellation makes sure opposite faces are divided along the same line
See figure http://www.ics.uci.edu/~eppstein/projects/tetra/
4, 5, 6, 0
0, 1, 5, 6
0, 2, 1, 6
4, 6, 7, 0
0, 7, 3, 6
0, 3, 2, 6
Introduces the additional diagonal edges in the Hex from 12 up to and including 18:
0 2 // 12 NEG_K
0 5 // 13 NEG_J
1 6 // 14 POS_I
3 6 // 15 POS_J
0 7 // 16 NEG_I
4 6 // 17 POS_K
0 6 // 18 Internal Diagonal
/// 4-----4------5
/// /| /| k POS_I = 0
/// 7 8 17 5 9 | NEG_I = 1
/// / | 13 / | | POS_J = 2
/// 7------6-----6 14| | NEG_J = 3
/// |16 0-----0--|---1 *------i POS_K = 4
/// 11 / 15 10 / / NEG_K = 5
/// | 3 12 | 1 / NO_FACE = 6
/// |/ |/ j
/// 3------2-----2
*/
//--------------------------------------------------------------------------------------------------
int HexGridIntersectionTools::planeHexIntersectionMCTet( const cvf::Plane& plane,
const cvf::Vec3d cell[8],
const size_t hexCornersIds[8],
std::vector<ClipVx>* triangleVxes,
std::vector<int>* cellFaceForEachTriangleEdge )
{
std::array<double, 8> cellCornerSqDistToPlane =
{
plane.distanceSquared( cell[0] ),
plane.distanceSquared( cell[1] ),
plane.distanceSquared( cell[2] ),
plane.distanceSquared( cell[3] ),
plane.distanceSquared( cell[4] ),
plane.distanceSquared( cell[5] ),
plane.distanceSquared( cell[6] ),
plane.distanceSquared( cell[7] ),
};
int cubeIndex = 0;
if (cellCornerSqDistToPlane[0] < 0) cubeIndex |= 1;
if (cellCornerSqDistToPlane[1] < 0) cubeIndex |= 2;
if (cellCornerSqDistToPlane[2] < 0) cubeIndex |= 4;
if (cellCornerSqDistToPlane[3] < 0) cubeIndex |= 8;
if (cellCornerSqDistToPlane[4] < 0) cubeIndex |= 16;
if (cellCornerSqDistToPlane[5] < 0) cubeIndex |= 32;
if (cellCornerSqDistToPlane[6] < 0) cubeIndex |= 64;
if (cellCornerSqDistToPlane[7] < 0) cubeIndex |= 128;
if (cubeIndex == 0 || cubeIndex == 255) return 0;
int tetCount = 0;
tetCount += planeMcTetIntersection(plane, cell, hexCornersIds, cellCornerSqDistToPlane.data(), { 4, 5, 6, 0 }, triangleVxes, cellFaceForEachTriangleEdge );
tetCount += planeMcTetIntersection(plane, cell, hexCornersIds, cellCornerSqDistToPlane.data(), { 0, 1, 5, 6 }, triangleVxes, cellFaceForEachTriangleEdge );
tetCount += planeMcTetIntersection(plane, cell, hexCornersIds, cellCornerSqDistToPlane.data(), { 0, 2, 1, 6 }, triangleVxes, cellFaceForEachTriangleEdge );
tetCount += planeMcTetIntersection(plane, cell, hexCornersIds, cellCornerSqDistToPlane.data(), { 4, 6, 7, 0 }, triangleVxes, cellFaceForEachTriangleEdge );
tetCount += planeMcTetIntersection(plane, cell, hexCornersIds, cellCornerSqDistToPlane.data(), { 0, 7, 3, 6 }, triangleVxes, cellFaceForEachTriangleEdge );
tetCount += planeMcTetIntersection(plane, cell, hexCornersIds, cellCornerSqDistToPlane.data(), { 0, 3, 2, 6 }, triangleVxes, cellFaceForEachTriangleEdge );
return tetCount;
}
//--------------------------------------------------------------------------------------------------
///
//--------------------------------------------------------------------------------------------------
cvf::uint HexGridIntersectionTools::planeMcTetIntersection( const cvf::Plane& plane,
const cvf::Vec3d hexCell[8],
const size_t hexCornersIds[8],
const double cornerDistToPlane[8],
const std::array<int, 4> & tetCell,
std::vector<ClipVx>* triangleVxes,
std::vector<int>* cellFaceForEachTriangleEdge )
{
static const int edgeEdgeCutsToCellFace[19][19] = {
// 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 4--------------4---------------5
{ 6, 5, 5, 5, 3, 6, 6, 6, 3, 3, 6, 6, 5, 3, 6, 6, 6, 6, 6 }, // 0 /|\__ __//|
{ 5, 6, 5, 5, 6, 0, 6, 6, 6, 0, 0, 6, 5, 6, 0, 6, 6, 6, 6 }, // 1 POS_I = 0 / | \___ _/ / |
{ 5, 5, 6, 5, 6, 6, 2, 6, 6, 6, 2, 2, 5, 6, 6, 2, 6, 6, 6 }, // 2 NEG_I = 1 7 | \__ __/ / | k
{ 5, 5, 5, 6, 6, 6, 6, 1, 1, 6, 6, 1, 5, 6, 6, 6, 1, 6, 6 }, // 3 POS_J = 2 / | 17___ _/ 5 | |
{ 3, 6, 6, 6, 6, 4, 4, 4, 3, 3, 6, 6, 6, 3, 6, 6, 6, 4, 6 }, // 4 NEG_J = 3 / 8 \___ / 9 |
{ 6, 0, 6, 6, 4, 6, 4, 4, 6, 0, 0, 6, 6, 6, 0, 6, 6, 4, 6 }, // 5 POS_K = 4 / | __/ \___ / | |
{ 6, 6, 2, 6, 4, 4, 6, 4, 6, 6, 2, 2, 6, 6, 6, 2, 6, 4, 6 }, // 6 NEG_K = 5 7---------------6--------------6 | *------i
{ 6, 6, 6, 1, 4, 4, 4, 6, 1, 6, 6, 1, 6, 6, 6, 6, 1, 4, 6 }, // 7 NO_FACE = 6 |\_ | __13 ____/_/|\_ | /
{ 3, 6, 6, 1, 3, 6, 6, 1, 6, 3, 6, 1, 6, 3, 6, 6, 1, 6, 6 }, // 8 | 16 | __/ _18_/ __/ | 14 | /
{ 3, 0, 6, 6, 3, 0, 6, 6, 3, 6, 0, 6, 6, 3, 0, 6, 6, 6, 6 }, // 9 | \_ | __/____/ _/ | \_ | j
{ 6, 0, 2, 6, 6, 0, 2, 6, 6, 0, 6, 2, 6, 6, 0, 2, 6, 6, 6 }, // 10 | \|__/__/ __/ | \|
{ 6, 6, 2, 1, 6, 6, 2, 1, 1, 6, 2, 6, 6, 6, 6, 2, 1, 6, 6 }, // 11 | 0-----------_/----0-----|------1
{ 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 }, // 12 11 / \__ __15 10 /
{ 3, 6, 6, 6, 3, 6, 6, 6, 3, 3, 6, 6, 6, 6, 6, 6, 6, 6, 6 }, // 13 | / \_ _/ | /
{ 6, 0, 6, 6, 6, 0, 6, 6, 6, 0, 0, 6, 6, 6, 6, 6, 6, 6, 6 }, // 14 | 3 __/ \___ | 1
{ 6, 6, 2, 6, 6, 6, 2, 6, 6, 6, 2, 2, 6, 6, 6, 6, 6, 6, 6 }, // 15 | / __/ 12__ | /
{ 6, 6, 6, 1, 6, 6, 6, 1, 1, 6, 6, 1, 6, 6, 6, 6, 6, 6, 6 }, // 16 | / __/ \___ | /
{ 6, 6, 6, 6, 4, 4, 4, 4, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 }, // 17 |/__/ \___|/
{ 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6 }, // 18 3---------------2--------------2
};
static const int cellCornerCellCornerToEdge[8][8] = {
// 0 1 2 3 4 5 6 7
{ -1, 0, 12, 3, 8, 13, 18, 16 }, // 0
{ 0, -1, 1, -1, -1, 9, 14, -1 }, // 1
{ 12, 1, -1, 2, -1, -1, 10, -1 }, // 2
{ 3, -1, 2, -1, -1, -1, 15, 11 }, // 3
{ 8, -1, -1, -1, -1, 4, 17, 7 }, // 4
{ 13, 9, -1, -1, 4, -1, 5, -1 }, // 5
{ 18, 14, 10, 15, 17, 5, -1, 6 }, // 6
{ 16, -1, -1, 11, 7, -1, 6, -1 }, // 7
};
cvf::uint ntri = 0;
int triindex = 0;
if( cornerDistToPlane[tetCell[0]] < 0 ) triindex |= 1;
if( cornerDistToPlane[tetCell[1]] < 0 ) triindex |= 2;
if( cornerDistToPlane[tetCell[2]] < 0 ) triindex |= 4;
if( cornerDistToPlane[tetCell[3]] < 0 ) triindex |= 8;
auto clipEdgeFunc = [&]( int hexCornerIdx0, int hexCornerIdx1 )
{
ClipVx cvx;
cvx.vx = planeLineIntersectionForMC( plane, hexCell[hexCornerIdx0], hexCell[hexCornerIdx1], &cvx.normDistFromEdgeVx1 );
cvx.clippedEdgeVx1Id = hexCornersIds[hexCornerIdx0];
cvx.clippedEdgeVx2Id = hexCornersIds[hexCornerIdx1];
return cvx;
};
auto addCellFaceStatusForTriangleEdges = [&]( int e11, int e12,
int e21, int e22,
int e31, int e32)
{
int cutEdge1 = cellCornerCellCornerToEdge[e11][e12];
int cutEdge2 = cellCornerCellCornerToEdge[e21][e22];
int cutEdge3 = cellCornerCellCornerToEdge[e31][e32];
CVF_ASSERT(cutEdge1 >= 0);
CVF_ASSERT(cutEdge2 >= 0);
CVF_ASSERT(cutEdge3 >= 0);
cellFaceForEachTriangleEdge->emplace_back( edgeEdgeCutsToCellFace[cutEdge1][cutEdge2] );
cellFaceForEachTriangleEdge->emplace_back( edgeEdgeCutsToCellFace[cutEdge2][cutEdge3] );
cellFaceForEachTriangleEdge->emplace_back( edgeEdgeCutsToCellFace[cutEdge3][cutEdge1] );
};
switch( triindex ) {
case 0x00:
case 0x0F:
break;
case 0x0E:
case 0x01:
{
triangleVxes->emplace_back( clipEdgeFunc( tetCell[0], tetCell[1] ) );
triangleVxes->emplace_back( clipEdgeFunc( tetCell[0], tetCell[2] ) );
triangleVxes->emplace_back( clipEdgeFunc( tetCell[0], tetCell[3] ) );
addCellFaceStatusForTriangleEdges(tetCell[0], tetCell[1],
tetCell[0], tetCell[2],
tetCell[0], tetCell[3]);
ntri++;
}
break;
case 0x0D:
case 0x02:
{
triangleVxes->emplace_back( clipEdgeFunc( tetCell[1], tetCell[0] ) );
triangleVxes->emplace_back( clipEdgeFunc( tetCell[1], tetCell[3] ) );
triangleVxes->emplace_back( clipEdgeFunc( tetCell[1], tetCell[2] ) );
addCellFaceStatusForTriangleEdges( tetCell[1], tetCell[0],
tetCell[1], tetCell[3],
tetCell[1], tetCell[2] );
ntri++;
}
break;
case 0x0C:
case 0x03:
{
triangleVxes->emplace_back( clipEdgeFunc( tetCell[0], tetCell[3] ) );
ClipVx vx1 = clipEdgeFunc( tetCell[0], tetCell[2] );
triangleVxes->push_back( vx1 );
ClipVx vx2 = clipEdgeFunc( tetCell[1], tetCell[3] );
triangleVxes->push_back( vx2 );
addCellFaceStatusForTriangleEdges( tetCell[0], tetCell[3],
tetCell[0], tetCell[2],
tetCell[1], tetCell[3] );
ntri++;
triangleVxes->push_back( vx2 );
triangleVxes->emplace_back( clipEdgeFunc( tetCell[1], tetCell[2] ) );
triangleVxes->push_back( vx1 );
addCellFaceStatusForTriangleEdges( tetCell[1], tetCell[3],
tetCell[1], tetCell[2],
tetCell[0], tetCell[2] );
ntri++;
}
break;
case 0x0B:
case 0x04:
{
triangleVxes->emplace_back( clipEdgeFunc( tetCell[2], tetCell[0] ) );
triangleVxes->emplace_back( clipEdgeFunc( tetCell[2], tetCell[1] ) );
triangleVxes->emplace_back( clipEdgeFunc( tetCell[2], tetCell[3] ) );
addCellFaceStatusForTriangleEdges( tetCell[2], tetCell[0],
tetCell[2], tetCell[1],
tetCell[2], tetCell[3] );
ntri++;
}
break;
case 0x0A:
case 0x05:
{
ClipVx vx0 = clipEdgeFunc( tetCell[0], tetCell[1] );
triangleVxes->push_back( vx0 );
ClipVx vx1 = clipEdgeFunc( tetCell[2], tetCell[3] );
triangleVxes->push_back( vx1 );
ClipVx vx2 = clipEdgeFunc( tetCell[0], tetCell[3] );
triangleVxes->push_back( vx2 );
addCellFaceStatusForTriangleEdges( tetCell[0], tetCell[1],
tetCell[2], tetCell[3],
tetCell[0], tetCell[3] );
ntri++;
triangleVxes->push_back( vx0 );
triangleVxes->emplace_back( clipEdgeFunc( tetCell[1], tetCell[2] ) );
triangleVxes->push_back( vx1 );
addCellFaceStatusForTriangleEdges( tetCell[0], tetCell[1],
tetCell[1], tetCell[2],
tetCell[2], tetCell[3] );
ntri++;
}
break;
case 0x09:
case 0x06:
{
ClipVx vx0 = clipEdgeFunc( tetCell[0], tetCell[1] );
triangleVxes->push_back( vx0 );
ClipVx vx1 = clipEdgeFunc( tetCell[1], tetCell[3] );
triangleVxes->push_back( vx1 );
ClipVx vx2 = clipEdgeFunc( tetCell[2], tetCell[3] );
triangleVxes->push_back( vx2 );
addCellFaceStatusForTriangleEdges( tetCell[0], tetCell[1],
tetCell[1], tetCell[3],
tetCell[2], tetCell[3] );
ntri++;
triangleVxes->push_back( vx0 );
triangleVxes->emplace_back( clipEdgeFunc( tetCell[0], tetCell[2] ) );
triangleVxes->push_back( vx2 );
addCellFaceStatusForTriangleEdges( tetCell[0], tetCell[1],
tetCell[0], tetCell[2],
tetCell[2], tetCell[3] );
ntri++;
}
break;
case 0x07:
case 0x08:
{
triangleVxes->emplace_back( clipEdgeFunc( tetCell[3], tetCell[0] ) );
triangleVxes->emplace_back( clipEdgeFunc( tetCell[3], tetCell[2] ) );
triangleVxes->emplace_back( clipEdgeFunc( tetCell[3], tetCell[1] ) );
addCellFaceStatusForTriangleEdges( tetCell[3], tetCell[0],
tetCell[3], tetCell[2],
tetCell[3], tetCell[1] );
ntri++;
}
break;
}
return ntri;
}
} // namespace cvf

View File

@ -92,6 +92,18 @@ public:
std::vector<ClipVx>* triangleVxes,
std::vector<int>* cellFaceForEachTriangleEdge);
static int planeHexIntersectionMCTet( const cvf::Plane& plane,
const cvf::Vec3d cell[8],
const size_t hexCornersIds[8],
std::vector<ClipVx>* triangleVxes,
std::vector<int>* cellFaceForEachTriangleEdge );
static cvf::uint planeMcTetIntersection( const cvf::Plane& plane,
const cvf::Vec3d hexCell[8],
const size_t hexCornersIds[8],
const double cornerDistToPlane[8],
const std::array<int, 4> & tetCell,
std::vector<ClipVx>* triangleVxes,
std::vector<int>* cellFaceForEachTriangleEdge );
};
}; // namespace caf