diff --git a/ApplicationCode/ModelVisualization/ModelVisualization_UnitTests/RivCellFaceIntersection-Test.cpp b/ApplicationCode/ModelVisualization/ModelVisualization_UnitTests/RivCellFaceIntersection-Test.cpp index 651688be9c..27dc2e840e 100644 --- a/ApplicationCode/ModelVisualization/ModelVisualization_UnitTests/RivCellFaceIntersection-Test.cpp +++ b/ApplicationCode/ModelVisualization/ModelVisualization_UnitTests/RivCellFaceIntersection-Test.cpp @@ -29,6 +29,7 @@ #include "cvfGeometryTools.h" #include "cvfBoundingBoxTree.h" +#include using namespace cvf; @@ -86,7 +87,7 @@ void ControlVolume::calculateCubeFaceStatus(const cvf::Vec3dArray& nodeCoords, d template -NodeType quadNormal (const ArrayWrapperConst& nodeCoords, +NodeType quadNormal (ArrayWrapperConst nodeCoords, const IndexType cubeFaceIndices[4] ) { return ( nodeCoords[cubeFaceIndices[2]] - nodeCoords[cubeFaceIndices[0]]) ^ @@ -143,7 +144,9 @@ public: private: QuadFaceIntersectorImplHandle * m_implementation; }; - +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- std::vector createVertices() { std::vector vxs; @@ -166,15 +169,221 @@ std::vector createVertices() return vxs; } - -std::vector > getCubeFaces() +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +std::vector getCubeFaces() { + std::vector cubeFaces; + + cvf::uint faces[4*4] = { + 0, 1, 2, 3, + 4, 5, 6, 7, + 5, 8, 9, 6, + 10, 11, 12, 13 + }; + + cubeFaces.resize(4); + cubeFaces[0] = &faces[0]; + cubeFaces[1] = &faces[4]; + cubeFaces[2] = &faces[8]; + cubeFaces[3] = &faces[12]; + + cubeFaces[3] = cubeFaces[2]; + return cubeFaces; +} + +std::ostream& operator<< (std::ostream& stream, std::vector v) +{ + for (size_t i = 0; i < v.size(); ++i) + { + stream << v[i] << " "; + } + return stream; +} + +//-------------------------------------------------------------------------------------------------- +/// +//-------------------------------------------------------------------------------------------------- +TEST(CellFaceIntersectionTst, Intersection1) +{ + std::vector nodes = createVertices(); + + std::vector additionalVertices; + + std::vector< std::vector > overlapPolygons; + std::vector faces = getCubeFaces(); + + EdgeIntersectStorage edgeIntersectionStorage; + edgeIntersectionStorage.setVertexCount(nodes.size()); + { + std::vector polygon; + bool isOk = false; + isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( + &polygon, + &additionalVertices, + edgeIntersectionStorage, + wrapArrayConst(&nodes), + faces[0].data(), + faces[1].data(), + 1e-6); + + EXPECT_EQ( 5, polygon.size()); + EXPECT_EQ( (size_t)2, additionalVertices.size()); + EXPECT_TRUE(isOk); + overlapPolygons.push_back(polygon); + std::cout << polygon << std::endl; + + } + + { + std::vector polygon; + bool isOk = false; + isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( + &polygon, + &additionalVertices, + edgeIntersectionStorage, + wrapArrayConst(&nodes), + faces[0].data(), + faces[2].data(), + 1e-6); + + EXPECT_EQ( 5, polygon.size()); + EXPECT_EQ( (size_t)4, additionalVertices.size()); + EXPECT_TRUE(isOk); + overlapPolygons.push_back(polygon); + std::cout << polygon << std::endl; + + } + + { + std::vector polygon; + bool isOk = false; + isOk = GeometryTools::calculateOverlapPolygonOfTwoQuads( + &polygon, + &additionalVertices, + edgeIntersectionStorage, + wrapArrayConst(&nodes), + faces[0].data(), + faces[3].data(), + 1e-6); + + EXPECT_EQ( 3, polygon.size()); + EXPECT_EQ( (size_t)6, additionalVertices.size()); + EXPECT_TRUE(isOk); + overlapPolygons.push_back(polygon); + std::cout << polygon << std::endl; + } + + nodes.insert(nodes.end(), additionalVertices.begin(), additionalVertices.end()); + std::vector basePolygon; + basePolygon.insert(basePolygon.begin(), faces[0].data(), &(faces[0].data()[4])); + + for (uint vxIdx = 0; vxIdx < nodes.size(); ++vxIdx) + { + bool inserted = GeometryTools::insertVertexInPolygon( + &basePolygon, + wrapArrayConst(&nodes), + vxIdx, + 1e-6 + ); + } + + EXPECT_EQ( 8, basePolygon.size()); + std::cout << "Bp: " << basePolygon << std::endl; + + for (size_t pIdx = 0; pIdx < overlapPolygons.size(); ++pIdx) + { + for (uint vxIdx = 0; vxIdx < nodes.size(); ++vxIdx) + { + bool inserted = GeometryTools::insertVertexInPolygon( + &overlapPolygons[pIdx], + wrapArrayConst(&nodes), + vxIdx, + 1e-6 + ); + } + + if (pIdx == 0) + { + EXPECT_EQ(5, overlapPolygons[pIdx].size()); + } + if (pIdx == 1) + { + EXPECT_EQ(5, overlapPolygons[pIdx].size()); + } + if (pIdx == 2) + { + EXPECT_EQ(4, overlapPolygons[pIdx].size()); + } + + std::cout << "Op" << pIdx << ":" << overlapPolygons[pIdx] << std::endl; + } + + + Vec3d normal = quadNormal(wrapArrayConst(&nodes), faces[0].data()); + std::vector faceOverlapPolygonWindingSameAsCubeFaceFlags; + faceOverlapPolygonWindingSameAsCubeFaceFlags.resize(overlapPolygons.size(), true); + + { + std::vector freeFacePolygon; + bool hasHoles = false; + + std::vector< std::vector* > overlapPolygonPtrs; + for (size_t pIdx = 0; pIdx < overlapPolygons.size(); ++pIdx) + { + overlapPolygonPtrs.push_back(&(overlapPolygons[pIdx])); + } + + GeometryTools::calculatePartiallyFreeCubeFacePolygon( + wrapArrayConst(&nodes), + wrapArrayConst(&basePolygon), + normal, + overlapPolygonPtrs, + faceOverlapPolygonWindingSameAsCubeFaceFlags, + &freeFacePolygon, + &hasHoles + ); + + EXPECT_EQ( 4, freeFacePolygon.size()); + EXPECT_FALSE(hasHoles); + std::cout << "FF1: " << freeFacePolygon << std::endl; + } + + { + std::vector freeFacePolygon; + bool hasHoles = false; + + std::vector< std::vector* > overlapPolygonPtrs; + for (size_t pIdx = 0; pIdx < 1; ++pIdx) + { + overlapPolygonPtrs.push_back(&(overlapPolygons[pIdx])); + } + + GeometryTools::calculatePartiallyFreeCubeFacePolygon( + wrapArrayConst(&nodes), + wrapArrayConst(&basePolygon), + normal, + overlapPolygonPtrs, + faceOverlapPolygonWindingSameAsCubeFaceFlags, + &freeFacePolygon, + &hasHoles + ); + + EXPECT_EQ( 9, freeFacePolygon.size()); + EXPECT_FALSE(hasHoles); + + std::cout << "FF2: " << freeFacePolygon << std::endl; + + } + } //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- + TEST(CellFaceIntersectionTst, Intersection) { diff --git a/ApplicationCode/ModelVisualization/cvfGeometryTools.cpp b/ApplicationCode/ModelVisualization/cvfGeometryTools.cpp index ba8a12a838..fbb10bae85 100644 --- a/ApplicationCode/ModelVisualization/cvfGeometryTools.cpp +++ b/ApplicationCode/ModelVisualization/cvfGeometryTools.cpp @@ -474,58 +474,6 @@ cvf::Vec3d GeometryTools::barycentricCoords(const cvf::Vec3d& t0, const cvf::Ve return m; } -//-------------------------------------------------------------------------------------------------- -/// Inserts the vertex into the polygon if it fits along one of the edges within the tolerance. -/// The method returns true if it was inserted, or if it was already in the polygon, or if it was -/// within the tolerance of an existing vertex in the polygon. -/// In the latter situation it replaces the previous vertex in the polygon. -/// -/// Todo: If a vertex is replaced, the VxToCv map in TimeStepGeometry should be updated -//-------------------------------------------------------------------------------------------------- -bool GeometryTools::insertVertexInPolygon(std::list >* polygon, const cvf::Vec3dArray& nodeCoords, cvf::uint vertexIndex, double tolerance) -{ - std::list >::iterator it; - for(it = polygon->begin(); it != polygon->end(); ++it) - { - if (it->first == vertexIndex) return true; - } - - -#if 1 - bool existsOrInserted = false; - for(it = polygon->begin(); it != polygon->end(); ++it) - { - if ( (nodeCoords[it->first] - nodeCoords[vertexIndex]).length() < tolerance) - { - if (vertexIndex < it->first) it->first = vertexIndex; - existsOrInserted = true; - } - } - - if (existsOrInserted) return true; -#endif - - // Insert vertex in polygon if the distance to one of the edges is small enough - - std::list >::iterator it2; - std::list >::iterator insertBefore; - - for (it = polygon->begin(); it != polygon->end(); ++it) - { - it2 = it; - ++it2; insertBefore = it2; if (it2 == polygon->end()) it2 = polygon->begin(); - - double sqDistToLine = GeometryTools::linePointSquareDist(nodeCoords[it->first], nodeCoords[it2->first], nodeCoords[vertexIndex]); - if (fabs(sqDistToLine) < tolerance*tolerance ) - { - it = polygon->insert(insertBefore, std::make_pair(vertexIndex, false)); - existsOrInserted = true; - } - } - - return existsOrInserted; -} - //-------------------------------------------------------------------------------------------------- /// //-------------------------------------------------------------------------------------------------- diff --git a/ApplicationCode/ModelVisualization/cvfGeometryTools.h b/ApplicationCode/ModelVisualization/cvfGeometryTools.h index 3951aeb6bd..ca57e7101a 100644 --- a/ApplicationCode/ModelVisualization/cvfGeometryTools.h +++ b/ApplicationCode/ModelVisualization/cvfGeometryTools.h @@ -37,9 +37,14 @@ public: LINES_OVERLAP }; - static bool insertVertexInPolygon(std::list >* polygon, const cvf::Vec3dArray& nodeCoords, cvf::uint vertexIndex, double tolerance); static void addMidEdgeNodes(std::list >* polygon, const cvf::Vec3dArray& nodes, EdgeSplitStorage& edgeSplitStorage, std::vector* createdVertexes); + template + static bool insertVertexInPolygon( std::vector * polygon, + ArrayWrapperConst nodeCoords, + IndexType vertexIndex, + double tolerance); + static IntersectionStatus inPlaneLineIntersect3D(const cvf::Vec3d& planeNormal, const cvf::Vec3d& p1, const cvf::Vec3d& p2, const cvf::Vec3d& p3, const cvf::Vec3d& p4, cvf::Vec3d* intersectionPoint, double* fractionAlongLine1, double* fractionAlongLine2, @@ -55,7 +60,8 @@ public: template - static bool calculateOverlapPolygonOfTwoQuads( std::vector * polygon, std::vector* createdVertexes, + static bool calculateOverlapPolygonOfTwoQuads( std::vector * polygon, + std::vector* createdVertexes, EdgeIntersectStorage& edgeIntersectionStorage, ArrayWrapperConst nodes, const IndexType cv1CubeFaceIndices[4], @@ -76,7 +82,7 @@ template class EdgeIntersectStorage { public: - void setVertexCount(IndexType size); + void setVertexCount(size_t size); bool findIntersection( IndexType e1P1, IndexType e1P2, IndexType e2P1, IndexType e2P2, IndexType* vxIndexIntersectionPoint, GeometryTools::IntersectionStatus* intersectionStatus, double* fractionAlongEdge1, double* fractionAlongEdge2); diff --git a/ApplicationCode/ModelVisualization/cvfGeometryTools.inl b/ApplicationCode/ModelVisualization/cvfGeometryTools.inl index 5f1d4615f1..c28b1f928b 100644 --- a/ApplicationCode/ModelVisualization/cvfGeometryTools.inl +++ b/ApplicationCode/ModelVisualization/cvfGeometryTools.inl @@ -3,6 +3,83 @@ namespace cvf { +//-------------------------------------------------------------------------------------------------- +/// Inserts the vertex into the polygon if it fits along one of the edges within the tolerance. +/// The method returns true if it was inserted, or if it was already in the polygon, or if it was +/// within the tolerance of an existing vertex in the polygon. +/// In the latter situation it replaces the previous vertex in the polygon. +/// +/// Todo: If a vertex is replaced, the VxToCv map in TimeStepGeometry should be updated +//-------------------------------------------------------------------------------------------------- +template +bool GeometryTools::insertVertexInPolygon( std::vector * polygon, + ArrayWrapperConst nodeCoords, + IndexType vertexIndex, + double tolerance) +{ + CVF_ASSERT(polygon); + + + // Check if vertex is directly included already + + for(std::vector::iterator it = polygon->begin(); it != polygon->end(); ++it) + { + if (*it == vertexIndex) return true; + } + +#if 1 + // Check if the new point is within tolerance of one of the polygon vertices + + bool existsOrInserted = false; + for(std::vector::iterator it = polygon->begin(); it != polygon->end(); ++it) + { + if ( (nodeCoords[*it] - nodeCoords[vertexIndex]).length() < tolerance) + { + if (vertexIndex < *it) *it = vertexIndex; + existsOrInserted = true; + } + } + + if (existsOrInserted) return true; +#endif + + + // Copy the start polygon to a list + + std::list listPolygon; + for (size_t pcIdx = 0; pcIdx < polygon->size(); ++pcIdx) + { + listPolygon.push_back((*polygon)[pcIdx]); + } + + // Insert vertex in polygon if the distance to one of the edges is small enough + + std::list::iterator it2; + std::list::iterator insertBefore; + + for (std::list::iterator it = listPolygon.begin(); it != listPolygon.end(); ++it) + { + it2 = it; + ++it2; insertBefore = it2; if (it2 == listPolygon.end()) it2 = listPolygon.begin(); + + double sqDistToLine = GeometryTools::linePointSquareDist(nodeCoords[*it], nodeCoords[*it2], nodeCoords[vertexIndex]); + if (fabs(sqDistToLine) < tolerance*tolerance ) + { + it = listPolygon.insert(insertBefore, vertexIndex); + existsOrInserted = true; + } + } + + // Write polygon back into the vector + + polygon->clear(); + for (std::list::iterator it = listPolygon.begin(); it != listPolygon.end(); ++it) + { + polygon->push_back(*it); + } + + return existsOrInserted; +} //-------------------------------------------------------------------------------------------------- @@ -116,6 +193,8 @@ bool GeometryTools::isPointTouchingIndexedPolygon( const cvf::Vec3d& polygonNor //-------------------------------------------------------------------------------------------------- /// Returns true if we get an actual polygon +/// The returned polygon will keep the winding from the first face. +/// The second face must have opposite winding of the first //-------------------------------------------------------------------------------------------------- template bool GeometryTools::calculateOverlapPolygonOfTwoQuads(std::vector * polygon, @@ -126,12 +205,12 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads(std::vector * p const IndexType cv2CubeFaceIndices[4], double tolerance) { + CVF_ASSERT(polygon); + CVF_ASSERT(createdVertexes); // Topology analysis - if (createdVertexes == NULL) return false; - - size_t newVertexIndex = nodes.size() + createdVertexes->size(); + IndexType newVertexIndex = static_cast(nodes.size() + createdVertexes->size()); bool cv1VxTouchCv2[4] = { false, false, false, false }; bool cv2VxTouchCv1[4] = { false, false, false, false }; @@ -387,13 +466,13 @@ bool GeometryTools::calculateOverlapPolygonOfTwoQuads(std::vector * p // Insert indices to the new intersection vertices into the polygon of // this connection according to fraction along edge - std::map sortingMap; + std::map sortingMap; for (i = 0; i < intersectionFractionsAlongEdge.size(); ++i) { sortingMap[intersectionFractionsAlongEdge[i]] = intersectionVxIndices[i]; } - std::map::iterator it; + std::map::iterator it; for (it = sortingMap.begin(); it != sortingMap.end(); ++it) { if (polygon->empty() || polygon->back() != it->second) @@ -478,7 +557,7 @@ void GeometryTools::calculatePartiallyFreeCubeFacePolygon(ArrayWrapperConst::const_iterator > VxIdxToPolygonPositionMap; + typedef std::map< IndexType, std::vector::const_iterator > VxIdxToPolygonPositionMap; CVF_ASSERT(m_partiallyFreeCubeFaceHasHoles); CVF_ASSERT(partialFacePolygon != NULL); @@ -745,7 +824,7 @@ void GeometryTools::calculatePartiallyFreeCubeFacePolygon(ArrayWrapperConst -void EdgeIntersectStorage::setVertexCount(IndexType size) +void EdgeIntersectStorage::setVertexCount(size_t size) { m_edgeIntsectMap.resize(size); } @@ -759,7 +838,7 @@ void EdgeIntersectStorage::canonizeAddress(IndexType& e1P1, IndexType flipE1E2 = (flipE1 ? e1P2: e1P1) > (flipE2 ? e2P2: e2P1); - static size_t temp; + static IndexType temp; if (flipE1) { temp = e1P1; @@ -835,15 +914,15 @@ bool EdgeIntersectStorage::findIntersection(IndexType e1P1, IndexType if (!m_edgeIntsectMap[e1P1].size()) return false; - std::map > >::iterator it; + std::map > >::iterator it; it = m_edgeIntsectMap[e1P1].find(e1P2); if (it == m_edgeIntsectMap[e1P1].end()) return false; - std::map >::iterator it2; + std::map >::iterator it2; it2 = it->second.find(e2P1); if (it2 == it->second.end()) return false; - std::map::iterator it3; + std::map::iterator it3; it3 = it2->second.find(e2P2); if (it3 == it2->second.end()) return false;