diff --git a/CMakePresets.json b/CMakePresets.json index 769d22e..3166b87 100644 --- a/CMakePresets.json +++ b/CMakePresets.json @@ -26,7 +26,6 @@ { "name": "gnu", "displayName": "GNU g++ compiler", - "generator": "Ninja", "inherits": "default" }, { diff --git a/src/core/Smoother.cpp b/src/core/Smoother.cpp index 9a15452..9dc38ea 100644 --- a/src/core/Smoother.cpp +++ b/src/core/Smoother.cpp @@ -8,6 +8,8 @@ #include #include +#include +#include #ifdef TESSELLATOR_EXECUTION_POLICIES #include @@ -19,6 +21,63 @@ namespace core { using namespace utils; using namespace meshTools; +namespace { + +using GridEntity = std::pair; + +struct OwnedCoordinates { + std::map edges; + std::map faces; + std::map interiors; +}; + +Cell canonicalCell(Cell cell, const SmootherTools& tools) +{ + for (Axis axis = X; axis <= Z; ++axis) { + cell[axis] = std::max(0, std::min(cell[axis], tools.numCellsDir(axis) - 1)); + } + return cell; +} + +OwnedCoordinates buildOwnedCoordinates( + const Coordinates& coordinates, + const SmootherTools& tools) +{ + OwnedCoordinates owned; + for (CoordinateId id = 0; id < coordinates.size(); ++id) { + const auto& coordinate = coordinates[id]; + if (tools.isRelativeInCellCorner(coordinate)) { + continue; // Grid corners are fixed throughout smoothing. + } + + Cell cell = canonicalCell(tools.toCell(coordinate), tools); + const auto edge = tools.getCellEdgeAxis(coordinate); + if (edge.first) { + owned.edges[{cell, edge.second}].insert(id); + continue; + } + const auto face = tools.getCellFaceAxis(coordinate); + if (face.first) { + owned.faces[{cell, face.second}].insert(id); + continue; + } + owned.interiors[cell].insert(id); + } + return owned; +} + +SmootherTools::IncidentElements buildIncidentElements(const Elements& elements) +{ + SmootherTools::IncidentElements incident; + for (const auto& element : elements) { + for (const auto id : element.vertices) { + incident[id].push_back(&element); + } + } + return incident; +} + +} Smoother::Smoother(const Mesh& mesh, const SmootherOptions& opts) : sT_(SmootherTools(mesh.grid)), @@ -35,6 +94,9 @@ Smoother::Smoother(const Mesh& mesh, const SmootherOptions& opts) : auto const singularIds = sT_.buildSingularIds(g.elements, mesh_.coordinates, opts_.featureDetectionAngle); + // Boundary remeshing is preprocessing. Keep its smooth sets separate + // from the ownership phases below: an element can touch several grid + // entities, while a coordinate has exactly one canonical owner. std::vector patchs; for (auto const& cell : sT_.buildCellElemMap(g.elements, mesh_.coordinates)) { for (auto const& p : @@ -46,34 +108,45 @@ Smoother::Smoother(const Mesh& mesh, const SmootherOptions& opts) : std::for_each(patchs.begin(), patchs.end(), [&](auto& p) { sT_.remeshBoundary(g.elements, res.coordinates, mesh_.coordinates, p); }); - - std::for_each(patchs.begin(), patchs.end(), [&](auto& p) { - sT_.collapsePointsOnCellEdges(res.coordinates, p, singularIds, opts_.contourAlignmentAngle); - }); - std::for_each( -#ifdef TESSELLATOR_EXECUTION_POLICIES - std::execution::par, -#endif - patchs.begin(), patchs.end(), [&](auto& p) { - sT_.collapsePointsOnCellFaces(res.coordinates, p, singularIds); - }); + const auto incidentElements = buildIncidentElements(g.elements); + const auto owned = buildOwnedCoordinates(res.coordinates, sT_); + IdSet edgeIds, faceIds, interiorIds; + for (const auto& entity : owned.edges) { + edgeIds.insert(entity.second.begin(), entity.second.end()); + } + for (const auto& entity : owned.faces) { + faceIds.insert(entity.second.begin(), entity.second.end()); + } + for (const auto& entity : owned.interiors) { + interiorIds.insert(entity.second.begin(), entity.second.end()); + } - std::for_each( -#ifdef TESSELLATOR_EXECUTION_POLICIES - std::execution::par, -#endif - patchs.begin(), patchs.end(), [&](auto& p) { - sT_.collapsePointsOnFeatureEdges(res.coordinates, p, singularIds); - }); + // Later phases may read an earlier boundary but never move it. + // This is intentionally sequential: entity ownership makes a future + // conflict graph possible, but entities sharing an element still need + // serialization until that scheduler is introduced. + for (const auto& patch : patchs) { + sT_.collapsePointsOnCellEdges(res.coordinates, patch, singularIds, + opts_.contourAlignmentAngle, edgeIds); + sT_.collapsePointsOnFeatureEdges(res.coordinates, patch, singularIds, + incidentElements, edgeIds); + } + meshTools::checkNoCellsAreCrossed(res); - std::for_each( -#ifdef TESSELLATOR_EXECUTION_POLICIES - std::execution::par, -#endif - patchs.begin(), patchs.end(), [&](auto& p) { - sT_.collapseInteriorPointsToBound(res.coordinates, p); - }); + for (const auto& patch : patchs) { + sT_.collapsePointsOnCellFaces(res.coordinates, patch, singularIds, faceIds); + sT_.collapsePointsOnFeatureEdges(res.coordinates, patch, singularIds, + incidentElements, faceIds); + } + meshTools::checkNoCellsAreCrossed(res); + + for (const auto& patch : patchs) { + sT_.collapsePointsOnFeatureEdges(res.coordinates, patch, singularIds, + incidentElements, interiorIds); + sT_.collapseInteriorPointsToBound(res.coordinates, patch, interiorIds); + } + meshTools::checkNoCellsAreCrossed(res); } @@ -95,4 +168,4 @@ Smoother::Smoother(const Mesh& mesh, const SmootherOptions& opts) : } } -} \ No newline at end of file +} diff --git a/src/core/SmootherTools.cpp b/src/core/SmootherTools.cpp index fd38246..dc46199 100644 --- a/src/core/SmootherTools.cpp +++ b/src/core/SmootherTools.cpp @@ -57,10 +57,66 @@ void SmootherTools::updateCoordinates( } } +bool SmootherTools::moveWouldCrossGrid( + const CoordinateId& id, + const Coordinate& destination, + const Coordinates& coordinates, + const IncidentElements& incidentElements) const +{ + const auto incident = incidentElements.find(id); + if (incident == incidentElements.end()) { + return false; + } + + for (const auto* element : incident->second) { + std::set commonCells; + bool firstVertex = true; + for (const auto vertexId : element->vertices) { + const auto& coordinate = vertexId == id ? destination : coordinates[vertexId]; + const auto touchingCells = getTouchingCells(coordinate); + if (firstVertex) { + commonCells = touchingCells; + firstVertex = false; + } + else { + for (auto cell = commonCells.begin(); cell != commonCells.end();) { + if (touchingCells.count(*cell) == 0) { + cell = commonCells.erase(cell); + } + else { + ++cell; + } + } + } + if (commonCells.empty()) { + return true; + } + } + } + return false; +} + +void SmootherTools::collapsePointsOnFeatureEdges( + Coordinates& coords, + const ElementsView& patch, + const SingularIds& singularIds, + const IdSet& movableIds) +{ + IncidentElements incidentElements; + for (const auto* element : patch) { + for (const auto id : element->vertices) { + incidentElements[id].push_back(element); + } + } + collapsePointsOnFeatureEdges(coords, patch, singularIds, incidentElements, movableIds); +} + void SmootherTools::collapsePointsOnFeatureEdges( Coordinates& coords, const ElementsView& patch, - const SingularIds& singularIds) + const SingularIds& singularIds, + const IncidentElements& incidentElements, + const IdSet& movableIds) { CoordGraph Point = CoordGraph(patch); CoordGraph edges = Point.getBoundaryGraph().intersect(singularIds.featureIds()); @@ -98,11 +154,31 @@ void SmootherTools::collapsePointsOnFeatureEdges( std::map toMove; for (auto const& i : validInterior) { + if (!movableIds.empty() && movableIds.count(i) == 0) { + continue; + } if (isRelativeInCellCorner(coords[i])) { continue; } - Coordinate closest = closestByDistance(coords, i, Point.getClosestVerticesInSet(i, validExterior)); + auto candidates = Point.getClosestVerticesInSet(i, validExterior); + const auto sourceCells = getTouchingCells(coords[i]); + for (auto candidate = candidates.begin(); candidate != candidates.end();) { + const auto candidateCells = getTouchingCells(coords[*candidate]); + const bool sharesCell = std::any_of( + candidateCells.begin(), candidateCells.end(), + [&](const Cell& cell) { return sourceCells.count(cell) != 0; }); + if (sharesCell) { + ++candidate; + } + else { + candidate = candidates.erase(candidate); + } + } + if (candidates.empty()) { + continue; + } + Coordinate closest = closestByDistance(coords, i, candidates); if (isRelativeInCellFace(coords[i]) && !areCoordOnSameFace(coords[i], closest)) { continue; } @@ -112,7 +188,12 @@ void SmootherTools::collapsePointsOnFeatureEdges( toMove[i] = closest; } - updateCoordinates(coords, toMove); + std::lock_guard lock(writingCoordinates_); + for (const auto& move : toMove) { + if (!moveWouldCrossGrid(move.first, move.second, coords, incidentElements)) { + coords[move.first] = move.second; + } + } } Coordinate SmootherTools::closestByDistance( @@ -169,7 +250,8 @@ void SmootherTools::collapsePointsOnCellEdges( Coordinates& coords, const ElementsView& patch, const SingularIds& singularIds, - double alignmentAngle) + double alignmentAngle, + const IdSet& movableIds) { { IdSet vertices = CoordGraph(patch).getVertices(); @@ -199,7 +281,10 @@ void SmootherTools::collapsePointsOnCellEdges( IdSet interiorValid = intersectWithIdSet(interior, protectedIds); - IdSet movable = classifyIds(interior, [&](auto i) {return !protectedIds.count(i); }).first; + IdSet movable = classifyIds(interior, [&](auto i) { + return !protectedIds.count(i) + && (movableIds.empty() || movableIds.count(i) != 0); + }).first; IdSet validIds = mergeIds(cG.getExterior(), interiorValid); std::map toMove; @@ -207,7 +292,11 @@ void SmootherTools::collapsePointsOnCellEdges( if (isRelativeInCellCorner(coords[i])) { continue; } - Coordinate closest = closestByDistance(coords, i, cG.getClosestVerticesInSet(i, validIds)); + const auto candidates = cG.getClosestVerticesInSet(i, validIds); + if (candidates.empty()) { + continue; + } + Coordinate closest = closestByDistance(coords, i, candidates); if (isRelativeInCellFace(coords[i]) && !areCoordOnSameFace(coords[i], closest)) { continue; } @@ -241,7 +330,8 @@ IdSet SmootherTools::getClosestValidByDistanceInCycle( void SmootherTools::collapsePointsOnCellFaces( Coordinates& coords, const ElementsView& patch, - const SingularIds& sIds) + const SingularIds& sIds, + const IdSet& movableIds) { std::map toMove; std::map > cyclesToValidOrOnFace; @@ -260,6 +350,9 @@ void SmootherTools::collapsePointsOnCellFaces( const IdSet& valid = kv.second.first; const IdSet& onCellFace = kv.second.second; for (auto const& id : onCellFace) { + if (!movableIds.empty() && movableIds.count(id) == 0) { + continue; + } try { IdSet candidates = getClosestValidByDistanceInCycle(id, cycle, valid); if (!candidates.empty()) { @@ -426,7 +519,7 @@ CoordGraph::Path SmootherTools::pathFromIdToAnyTarget( } const std::size_t i = it - cycle.begin(); - CoordGraph::Path res(startId); + CoordGraph::Path res{startId}; for (std::size_t d = 1; d < cycle.size(); d++) { CoordinateId idTest; if (forward) { @@ -485,13 +578,17 @@ Coordinates SmootherTools::collapsePointsOnContour( void SmootherTools::collapseInteriorPointsToBound( Coordinates& coords, - const ElementsView& patch) + const ElementsView& patch, + const IdSet& movableIds) { IdSet bound, interior; std::tie(bound, interior) = CoordGraph(patch).getBoundAndInteriorVertices(); std::map toMove; for (auto const& vI : interior) { + if (!movableIds.empty() && movableIds.count(vI) == 0) { + continue; + } toMove[vI] = closestByDistance(coords, vI, bound); } @@ -573,4 +670,4 @@ void SmootherTools::reorientSingleElement( } } -} \ No newline at end of file +} diff --git a/src/core/SmootherTools.h b/src/core/SmootherTools.h index 973de4f..3ab957b 100644 --- a/src/core/SmootherTools.h +++ b/src/core/SmootherTools.h @@ -13,6 +13,8 @@ namespace core { class SmootherTools : public utils::GridTools { public: + using IncidentElements = std::map; + class SingularIds { public: SingularIds(const IdSet& featureIds, const IdSet& contourIds, const IdSet& cornerIds) : @@ -39,7 +41,15 @@ class SmootherTools : public utils::GridTools { void collapsePointsOnFeatureEdges( Coordinates& res, const ElementsView& patch, - const SingularIds& singularIds); + const SingularIds& singularIds, + const IdSet& movableIds = {}); + + void collapsePointsOnFeatureEdges( + Coordinates& res, + const ElementsView& patch, + const SingularIds& singularIds, + const IncidentElements& incidentElements, + const IdSet& movableIds = {}); Coordinates collapsePointsOnContour( const Elements& elems, @@ -50,12 +60,14 @@ class SmootherTools : public utils::GridTools { Coordinates& res, const ElementsView& patch, const SingularIds& singularIds , - double alignmentAngle); + double alignmentAngle, + const IdSet& movableIds = {}); void collapsePointsOnCellFaces( Coordinates& res, const ElementsView& patch, - const SingularIds&); + const SingularIds&, + const IdSet& movableIds = {}); void remeshBoundary( Elements& es, @@ -70,7 +82,8 @@ class SmootherTools : public utils::GridTools { void collapseInteriorPointsToBound( Coordinates& coords, - const ElementsView& patch); + const ElementsView& patch, + const IdSet& movableIds = {}); void remeshElementsToOneInteriorPoint( Elements& es, @@ -83,11 +96,19 @@ class SmootherTools : public utils::GridTools { const ElementsView& patch); private: + friend class SmootherToolsTestAccess; + std::mutex writingCoordinates_; std::mutex writingElements_; void updateCoordinates(Coordinates& res, std::map toMove); + bool moveWouldCrossGrid( + const CoordinateId& id, + const Coordinate& destination, + const Coordinates& coordinates, + const IncidentElements& incidentElements) const; + static CoordinateId getClosestEndOfPaths( const std::vector& paths); @@ -136,4 +157,4 @@ class SmootherTools : public utils::GridTools { }; } -} \ No newline at end of file +} diff --git a/src/utils/CoordGraph.cpp b/src/utils/CoordGraph.cpp index e4242f6..0e5e0c3 100644 --- a/src/utils/CoordGraph.cpp +++ b/src/utils/CoordGraph.cpp @@ -81,6 +81,10 @@ CoordGraph::CoordGraph(const Elements& elems) CoordGraph::CoordGraph(const ElementsView& es) { for (auto const& e : es) { + if (e->vertices.size() == 1) { + addVertex(e->vertices.front()); + continue; + } for (std::size_t i = 0; i < e->vertices.size(); i++) { this->addEdge( e->vertices[i], @@ -95,6 +99,10 @@ CoordGraph::CoordGraph(const ElementsView& es) CoordGraph::CoordGraph(const Paths& paths) { for (const auto& p : paths) { + if (p.size() == 1) { + addVertex(p.front()); + continue; + } for (std::size_t i = 0; i < p.size(); i++) { this->addEdge( p[i], @@ -221,7 +229,7 @@ IdSet CoordGraph::getClosestVerticesInSet( } auto path = findShortestPath(vI, vB); if (path.empty()) { - throw std::runtime_error("Can not find path to point in set."); + continue; } paths.insert(path); } @@ -586,4 +594,4 @@ std::vector CoordGraph::findCycles() const } } -} \ No newline at end of file +} diff --git a/test/core/SmootherToolsTest.cpp b/test/core/SmootherToolsTest.cpp index 78e2670..e93a9d0 100644 --- a/test/core/SmootherToolsTest.cpp +++ b/test/core/SmootherToolsTest.cpp @@ -419,6 +419,30 @@ class SmootherToolsTest : public ::testing::Test { const double alignmentAngle = 5.0; }; +class SmootherToolsTestAccess { +public: + static CoordGraph::Path pathFromIdToAnyTarget( + const CoordinateId startId, + const CoordGraph::Path& cycle, + const bool forward, + const IdSet& target) + { + return SmootherTools::pathFromIdToAnyTarget(startId, cycle, forward, target); + } +}; + +TEST_F(SmootherToolsTest, pathFromLargeCoordinateIdContainsTheIdOnce) +{ + const CoordinateId startId = 1000; + const CoordinateId targetId = 1001; + const CoordGraph::Path cycle = {999, startId, targetId}; + + const auto path = SmootherToolsTestAccess::pathFromIdToAnyTarget( + startId, cycle, true, {targetId}); + + EXPECT_EQ((CoordGraph::Path{startId, targetId}), path); +} + /// Remesh patch with interior points to remove them. /// \verbatim /// 1 ---- 2 @@ -673,6 +697,105 @@ TEST_F(SmootherToolsTest, collapsePointsOnFeatureEdges_singlePatch) EXPECT_EQ(7, countDifferentCoordinates(cs)); } +TEST_F(SmootherToolsTest, collapsePointsOnFeatureEdges_withDisconnectedTargets) +{ + Mesh m; + m.grid = buildUnitLengthGrid(1.0); + m.coordinates = { + Coordinate({0.1, 0.1, 0.2}), + Coordinate({0.5, 0.1, 0.2}), + Coordinate({0.9, 0.1, 0.2}), + Coordinate({0.5, 0.4, 0.2}), + Coordinate({0.1, 0.6, 0.8}), + Coordinate({0.5, 0.6, 0.8}), + Coordinate({0.9, 0.6, 0.8}), + Coordinate({0.5, 0.9, 0.8}), + }; + m.groups = {Group()}; + m.groups[0].elements = { + Element({0, 1, 3}), + Element({1, 2, 3}), + Element({4, 5, 7}), + Element({5, 6, 7}), + }; + + Elements& elements = m.groups[0].elements; + const ElementsView disconnectedPatch = { + &elements[0], &elements[1], &elements[2], &elements[3] + }; + const SmootherTools::SingularIds singularIds( + {0, 1, 2, 4, 5, 6}, {}, {}); + + EXPECT_NO_THROW(SmootherTools(m.grid).collapsePointsOnFeatureEdges( + m.coordinates, disconnectedPatch, singularIds)); + EXPECT_EQ(m.coordinates[0], m.coordinates[1]); + EXPECT_EQ(m.coordinates[4], m.coordinates[5]); +} + +TEST_F(SmootherToolsTest, collapsePointsOnFeatureEdgesDoesNotCrossTouchingCell) +{ + Mesh m; + m.grid = buildUnitLengthGrid(0.25); + m.coordinates = { + Coordinate({0.9, 0.1, 0.2}), + Coordinate({1.2, 0.1, 0.2}), + Coordinate({1.9, 0.1, 0.2}), + Coordinate({1.5, 0.4, 0.2}), + }; + m.groups = {Group()}; + m.groups[0].elements = { + Element({0, 1, 3}), + Element({1, 2, 3}), + }; + + Elements& elements = m.groups[0].elements; + const ElementsView patch = {&elements[0], &elements[1]}; + const SmootherTools::SingularIds singularIds({0, 1, 2}, {}, {}); + + SmootherTools(m.grid).collapsePointsOnFeatureEdges( + m.coordinates, patch, singularIds); + + EXPECT_EQ(Coordinate({1.2, 0.1, 0.2}), m.coordinates[1]); + EXPECT_FALSE(GridTools(m.grid).elementCrossesGrid(elements[1], m.coordinates)); +} + +TEST_F(SmootherToolsTest, collapsePointsOnFeatureEdgesRejectsCombinedMovesThatCrossGrid) +{ + Mesh m; + m.grid = buildUnitLengthGrid(0.25); + m.coordinates = { + Coordinate({1.0, 0.2, 0.2}), + Coordinate({1.0, 1.0, 0.2}), + Coordinate({1.0, 1.0, 0.8}), + Coordinate({1.0, 1.8, 0.8}), + Coordinate({0.8, 1.0, 0.5}), + }; + m.groups = {Group()}; + m.groups[0].elements = { + Element({0, 1, 4}), + Element({1, 2, 4}), + Element({2, 3, 4}), + }; + + Elements& elements = m.groups[0].elements; + const ElementsView patch = {&elements[0], &elements[1], &elements[2]}; + const SmootherTools::SingularIds singularIds({0, 1, 2, 3}, {}, {}); + SmootherTools::IncidentElements incidentElements; + for (const auto& element : elements) { + for (const auto id : element.vertices) { + incidentElements[id].push_back(&element); + } + } + + SmootherTools(m.grid).collapsePointsOnFeatureEdges( + m.coordinates, patch, singularIds, incidentElements); + + EXPECT_EQ(m.coordinates[0], m.coordinates[1]); + EXPECT_NE(m.coordinates[3], m.coordinates[2]); + EXPECT_FALSE(GridTools(m.grid).elementCrossesGrid( + elements[1], m.coordinates)); +} + TEST_F(SmootherToolsTest, collapsePointsOnFeatureEdges_threePatches) { Mesh m = buildCornerMesh(); @@ -860,4 +983,4 @@ TEST_F(SmootherToolsTest, collapsePointsOnFeatureEdges_feature_in_interior) EXPECT_EQ(collapsed[5], collapsed[4]); } -} \ No newline at end of file +} diff --git a/test/utils/CoordGraphTest.cpp b/test/utils/CoordGraphTest.cpp index 6a98958..54e9259 100644 --- a/test/utils/CoordGraphTest.cpp +++ b/test/utils/CoordGraphTest.cpp @@ -341,6 +341,20 @@ TEST_F(CoordGraphTest, ctors) EXPECT_EQ(cG.verticesSize(), cGViews.verticesSize()); EXPECT_EQ(cG.edgesSize(), cGViews.edgesSize()); } + +TEST_F(CoordGraphTest, constructorsAcceptSingleVertexInputs) +{ + const Elements elements = { + Element({7}, Element::Type::Node), + }; + const CoordGraph fromElements(elements); + const CoordGraph fromPaths(CoordGraph::Paths{{9}}); + + EXPECT_EQ(IdSet({7}), fromElements.getVertices()); + EXPECT_EQ(IdSet({9}), fromPaths.getVertices()); + EXPECT_EQ(0, fromElements.edgesSize()); + EXPECT_EQ(0, fromPaths.edgesSize()); +} TEST_F(CoordGraphTest, adding_edge_with_non_existingvertices) { CoordGraph g; @@ -563,6 +577,17 @@ TEST_F(CoordGraphTest, getClosestVerticesInSet_2) EXPECT_EQ(IdSet({ 3 }), g.getClosestVerticesInSet(2, { 3, 4 })); } + +TEST_F(CoordGraphTest, getClosestVerticesInSet_ignoresDisconnectedTargets) +{ + CoordGraph g; + g.addEdge(1, 2); + g.addEdge(3, 4); + + EXPECT_EQ(IdSet({2}), g.getClosestVerticesInSet(1, {2, 3})); + EXPECT_TRUE(g.getClosestVerticesInSet(1, {3}).empty()); +} + TEST_F(CoordGraphTest, graphIntersection) { @@ -1341,4 +1366,4 @@ TEST_F(CoordGraphTest, difference) } } -} \ No newline at end of file +}