Fix/conformal smoothing

CMakePresets.json

@@ -26,7 +26,6 @@
         {
             "name": "gnu",
             "displayName": "GNU g++ compiler",
-            "generator": "Ninja",
             "inherits": "default"
         },
         {

src/core/Smoother.cpp

@@ -8,6 +8,8 @@
 
 #include <assert.h>
 #include <algorithm>
+#include <map>
+#include <set>
 
 #ifdef TESSELLATOR_EXECUTION_POLICIES
 #include <execution>
@@ -19,6 +21,63 @@ namespace core {
 using namespace utils;
 using namespace meshTools;
 
+namespace {
+
+using GridEntity = std::pair<Cell, Axis>;
+
+struct OwnedCoordinates {
+    std::map<GridEntity, IdSet> edges;
+    std::map<GridEntity, IdSet> faces;
+    std::map<Cell, IdSet> interiors;
+};
+
+Cell canonicalCell(Cell cell, const SmootherTools& tools)
+{
+    for (Axis axis = X; axis <= Z; ++axis) {
+        cell[axis] = std::max<CellDir>(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<ElementsView> 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) :
 }
 
 }
-}
+}

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<Cell> 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<CoordinateId, Coordinate> 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<std::mutex> 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,15 +281,22 @@ 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<CoordinateId, Coordinate> toMove;
         for (auto const& i : movable) {
             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<CoordinateId, Coordinate> toMove;
     std::map <CoordGraph::Path, std::pair<IdSet, IdSet>> 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<CoordinateId, Coordinate> 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(
 }
 
 }
-}
+}