From 37f684386f635b44973afcdb640bc7a72138afd5 Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Tue, 21 Oct 2025 16:00:05 +0200 Subject: [PATCH 01/23] Adds json support for conformal mesher in launcher. Modifies conformalMesher to add selective structuring of non-conformal cells after meshing --- src/app/launcher.cpp | 52 +++++++++++++++++++++++++++++++-- src/app/launcher.h | 4 +++ src/meshers/ConformalMesher.cpp | 14 +++++++-- test/app/launcherTest.cpp | 39 ++++++++++++++++++++++++- 4 files changed, 103 insertions(+), 6 deletions(-) diff --git a/src/app/launcher.cpp b/src/app/launcher.cpp index 42e2805..596b223 100644 --- a/src/app/launcher.cpp +++ b/src/app/launcher.cpp @@ -2,6 +2,7 @@ #include "vtkIO.h" #include "meshers/StructuredMesher.h" +#include "meshers/ConformalMesher.h" #include "utils/GridTools.h" #include @@ -11,12 +12,13 @@ #include #include #include - +#include namespace meshlib::app { using namespace vtkIO; + namespace po = boost::program_options; Grid parseGridFromJSON(const nlohmann::json &j) @@ -61,6 +63,49 @@ Mesh readMesh(const std::string &fn) return res; } +std::string readMesherType(const std::string &fn) +{ + nlohmann::json j; + { + std::ifstream i(fn); + i >> j; + } + if (j["mesher"].contains("type")) { + return j["mesher"]["type"]; + } else { + return meshlib::app::structured_mesher; + } +} + +meshlib::meshers::ConformalMesherOptions readConformalMesherOptions(const std::string &fn) +{ + nlohmann::json j; + { + std::ifstream i(fn); + i >> j; + } + meshlib::meshers::ConformalMesherOptions res; + if (j["mesher"].contains("options")) { + res.snapperOptions.edgePoints = j["mesher"]["options"]["edgePoints"]; + res.snapperOptions.forbiddenLength = j["mesher"]["options"]["forbiddenLength"]; + } else { + res.snapperOptions.edgePoints = 4; + res.snapperOptions.forbiddenLength = 0; + } + return res; +} +std::unique_ptr buildMesher(const Mesh &in, const std::string &fn) +{ + auto mesherType = readMesherType(fn); + if (mesherType == meshlib::app::structured_mesher) { + return std::make_unique(meshlib::meshers::StructuredMesher{in}); + } else if (mesherType == meshlib::app::conformal_mesher) { + return std::make_unique(meshlib::meshers::ConformalMesher{in, readConformalMesherOptions(fn)}); + } else { + throw std::runtime_error("Unsupported mesher type"); + } +} + int launcher(int argc, const char* argv[]) { po::options_description desc("Allowed options"); @@ -84,9 +129,10 @@ int launcher(int argc, const char* argv[]) Mesh mesh = readMesh(inputFilename); + // Mesh - meshlib::meshers::StructuredMesher mesher{mesh}; - Mesh resultMesh = mesher.mesh(); + auto mesher = buildMesher(mesh, inputFilename); + Mesh resultMesh = mesher->mesh(); std::filesystem::path outputFolder = getFolder(inputFilename); auto basename = getBasename(inputFilename); diff --git a/src/app/launcher.h b/src/app/launcher.h index 4e33199..ca5388b 100644 --- a/src/app/launcher.h +++ b/src/app/launcher.h @@ -1,9 +1,13 @@ #pragma once #include "types/Mesh.h" +#include namespace meshlib::app { +const std::string conformal_mesher ("conformal"); +const std::string structured_mesher ("structured"); + int launcher(int argc, const char* argv[]); } \ No newline at end of file diff --git a/src/meshers/ConformalMesher.cpp b/src/meshers/ConformalMesher.cpp index 6811c55..5406679 100644 --- a/src/meshers/ConformalMesher.cpp +++ b/src/meshers/ConformalMesher.cpp @@ -5,10 +5,12 @@ #include "core/Collapser.h" #include "core/Smoother.h" #include "core/Snapper.h" +#include "core/Staircaser.h" #include "utils/GridTools.h" #include "utils/MeshTools.h" #include "utils/CoordGraph.h" +#include "utils/RedundancyCleaner.h" namespace meshlib::meshers { @@ -193,10 +195,18 @@ Mesh ConformalMesher::mesh() const log("Non-conformal cells found: " + std::to_string(nonConformalCells.size()), 1); // Calls structurer to mesh only those cells. - + log("Structuring non-conformal cells.", 1); + res = Staircaser{ res }.getSelectiveMesh(nonConformalCells,Staircaser::GapsFillingType::Split); + RedundancyCleaner::removeOverlappedDimensionOneAndLowerElementsAndEquivalentSurfaces(res); + logNumberOfTriangles(countMeshElementsIf(res, isTriangle)); + + // Find cells which break conformal FDTD rules after selective structuring. + nonConformalCells = findNonConformalCells(res); + log("Non-conformal cells found after Selective Structuring: " + std::to_string(nonConformalCells.size()), 1); + + // Merges triangles which are on same cell face. - reduceGrid(res, originalGrid_); // Converts relatives to absolutes. diff --git a/test/app/launcherTest.cpp b/test/app/launcherTest.cpp index 4188f2d..bcbd5c2 100644 --- a/test/app/launcherTest.cpp +++ b/test/app/launcherTest.cpp @@ -24,6 +24,15 @@ TEST_F(LauncherTest, launches_alhambra_case) EXPECT_EQ(exitCode, EXIT_SUCCESS); } +TEST_F(LauncherTest, launches_conformal_alhambra_case) +{ + int ac = 3; + const char* av[] = { NULL, "-i", "testData/cases/alhambra/alhambra.conformal.tessellator.json"}; + int exitCode; + EXPECT_NO_THROW(exitCode = meshlib::app::launcher(ac, av)); + EXPECT_EQ(exitCode, EXIT_SUCCESS); +} + TEST_F(LauncherTest, launches_sphere_case) { @@ -34,6 +43,24 @@ TEST_F(LauncherTest, launches_sphere_case) EXPECT_EQ(exitCode, EXIT_SUCCESS); } +TEST_F(LauncherTest, launches_conformal_sphere_case) +{ + int ac = 3; + const char* av[] = { NULL, "-i", "testData/cases/sphere/sphere.conformal.tessellator.json"}; + int exitCode; + EXPECT_NO_THROW(exitCode = meshlib::app::launcher(ac, av)); + EXPECT_EQ(exitCode, EXIT_SUCCESS); +} + +TEST_F(LauncherTest, launches_conformal_thinCylinder_case) +{ + int ac = 3; + const char* av[] = { NULL, "-i", "testData/cases/thinCylinder/thinCylinder.conformal.tessellator.json"}; + int exitCode; + EXPECT_NO_THROW(exitCode = meshlib::app::launcher(ac, av)); + EXPECT_EQ(exitCode, EXIT_SUCCESS); +} + TEST_F(LauncherTest, launches_thinCylinder_case) { int ac = 3; @@ -50,4 +77,14 @@ TEST_F(LauncherTest, launches_cone_case) int exitCode; EXPECT_NO_THROW(exitCode = meshlib::app::launcher(ac, av)); EXPECT_EQ(exitCode, EXIT_SUCCESS); -} \ No newline at end of file +} + +TEST_F(LauncherTest, launches_conformal_cone_case) +{ + int ac = 3; + const char* av[] = { NULL, "-i", "testData/cases/cone/cone.conformal.tessellator.json"}; + int exitCode; + EXPECT_NO_THROW(exitCode = meshlib::app::launcher(ac, av)); + EXPECT_EQ(exitCode, EXIT_SUCCESS); +} + From 35377db8da30c3bc49d007a71f37373bcd130471 Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Wed, 22 Oct 2025 08:28:05 +0200 Subject: [PATCH 02/23] If conformalMesher options are not given, default are taken from the SnapperOptions struct --- src/app/launcher.cpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/app/launcher.cpp b/src/app/launcher.cpp index 596b223..0692747 100644 --- a/src/app/launcher.cpp +++ b/src/app/launcher.cpp @@ -88,9 +88,6 @@ meshlib::meshers::ConformalMesherOptions readConformalMesherOptions(const std::s if (j["mesher"].contains("options")) { res.snapperOptions.edgePoints = j["mesher"]["options"]["edgePoints"]; res.snapperOptions.forbiddenLength = j["mesher"]["options"]["forbiddenLength"]; - } else { - res.snapperOptions.edgePoints = 4; - res.snapperOptions.forbiddenLength = 0; } return res; } From 45d99b32a5148e9aeb67090f8520309c2caa180d Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Wed, 22 Oct 2025 08:33:28 +0200 Subject: [PATCH 03/23] Adds conformal json files --- .../alhambra.conformal.tessellator.json | 17 +++++++++++++++++ .../cases/cone/cone.conformal.tessellator.json | 17 +++++++++++++++++ .../sphere/sphere.conformal.tessellator.json | 17 +++++++++++++++++ .../thinCylinder.conformal.tessellator.json | 17 +++++++++++++++++ 4 files changed, 68 insertions(+) create mode 100644 testData/cases/alhambra/alhambra.conformal.tessellator.json create mode 100644 testData/cases/cone/cone.conformal.tessellator.json create mode 100644 testData/cases/sphere/sphere.conformal.tessellator.json create mode 100644 testData/cases/thinCylinder/thinCylinder.conformal.tessellator.json diff --git a/testData/cases/alhambra/alhambra.conformal.tessellator.json b/testData/cases/alhambra/alhambra.conformal.tessellator.json new file mode 100644 index 0000000..f2b6c20 --- /dev/null +++ b/testData/cases/alhambra/alhambra.conformal.tessellator.json @@ -0,0 +1,17 @@ +{ + "grid": { + "numberOfCells": [60, 60, 10], + "boundingBox": [ + [-60.0, -60.0, -10.0], + [ 60.0, 60.0, 10.0] + ] + }, + "object": {"filename": "alhambra.stl"}, + "mesher": { + "type" : "conformal", + "options" : { + "edgePoints" : 6, + "forbiddenLength" : 0.15 + } + } +} \ No newline at end of file diff --git a/testData/cases/cone/cone.conformal.tessellator.json b/testData/cases/cone/cone.conformal.tessellator.json new file mode 100644 index 0000000..dac46e0 --- /dev/null +++ b/testData/cases/cone/cone.conformal.tessellator.json @@ -0,0 +1,17 @@ +{ + "grid": { + "numberOfCells": [40, 40, 120], + "boundingBox": [ + [-2, -2, -1], + [ 2, 2, 11] + ] + }, + "object": {"filename": "cone.stl"}, + "mesher": { + "type" : "conformal", + "options" : { + "edgePoints" : 4, + "forbiddenLength" : 0.25 + } + } +} \ No newline at end of file diff --git a/testData/cases/sphere/sphere.conformal.tessellator.json b/testData/cases/sphere/sphere.conformal.tessellator.json new file mode 100644 index 0000000..67d7ef2 --- /dev/null +++ b/testData/cases/sphere/sphere.conformal.tessellator.json @@ -0,0 +1,17 @@ +{ + "grid": { + "numberOfCells": [50, 50, 50], + "boundingBox": [ + [-100.0, -100.0, -100.0], + [ 100.0, 100.0, 100.0] + ] + }, + "object": {"filename": "sphere.stl"}, + "mesher": { + "type" : "conformal", + "options" : { + "edgePoints" : 4, + "forbiddenLength" : 0.25 + } + } +} \ No newline at end of file diff --git a/testData/cases/thinCylinder/thinCylinder.conformal.tessellator.json b/testData/cases/thinCylinder/thinCylinder.conformal.tessellator.json new file mode 100644 index 0000000..d76cd2d --- /dev/null +++ b/testData/cases/thinCylinder/thinCylinder.conformal.tessellator.json @@ -0,0 +1,17 @@ +{ + "grid": { + "numberOfCells": [20, 20, 30], + "boundingBox": [ + [-1, -1, -1], + [ 1, 1, 2] + ] + }, + "object": {"filename": "thinCylinder.stl"}, + "mesher": { + "type" : "conformal", + "options" : { + "edgePoints" : 4, + "forbiddenLength" : 0.25 + } + } +} \ No newline at end of file From fe1a5864557f14b6a1e42e2da0be1fd592688bd2 Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Thu, 23 Oct 2025 18:30:58 +0200 Subject: [PATCH 04/23] Adds dynamic extension to the mesh file, depending on mesher used --- src/app/launcher.cpp | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/src/app/launcher.cpp b/src/app/launcher.cpp index 0692747..0c937af 100644 --- a/src/app/launcher.cpp +++ b/src/app/launcher.cpp @@ -63,6 +63,7 @@ Mesh readMesh(const std::string &fn) return res; } + std::string readMesherType(const std::string &fn) { nlohmann::json j; @@ -77,6 +78,18 @@ std::string readMesherType(const std::string &fn) } } +std::string readExtension(const std::string &fn) +{ + auto mesherType = readMesherType(fn); + if (mesherType == meshlib::app::structured_mesher) { + return "str"; + } else if (mesherType == meshlib::app::conformal_mesher) { + return "cmsh"; + } else { + throw std::runtime_error("Unsupported mesher type"); + } +} + meshlib::meshers::ConformalMesherOptions readConformalMesherOptions(const std::string &fn) { nlohmann::json j; @@ -133,7 +146,9 @@ int launcher(int argc, const char* argv[]) std::filesystem::path outputFolder = getFolder(inputFilename); auto basename = getBasename(inputFilename); - exportMeshToVTU(outputFolder / (basename + ".tessellator.str.vtk"), resultMesh); + auto extension = readExtension(inputFilename); + + exportMeshToVTU(outputFolder / (basename + ".tessellator." + extension + ".vtk"), resultMesh); exportGridToVTU(outputFolder / (basename + ".tessellator.grid.vtk"), resultMesh.grid); return EXIT_SUCCESS; From cb55b5c6c7afd6416e616f8b3db78545c5ec4309 Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Tue, 28 Oct 2025 09:20:24 +0100 Subject: [PATCH 05/23] Repo updates --- src/meshers/ConformalMesher.cpp | 2 +- tmp/sphere.tessellator.json | 10 ++++++++++ 2 files changed, 11 insertions(+), 1 deletion(-) create mode 100644 tmp/sphere.tessellator.json diff --git a/src/meshers/ConformalMesher.cpp b/src/meshers/ConformalMesher.cpp index 5406679..dca455c 100644 --- a/src/meshers/ConformalMesher.cpp +++ b/src/meshers/ConformalMesher.cpp @@ -196,7 +196,7 @@ Mesh ConformalMesher::mesh() const // Calls structurer to mesh only those cells. log("Structuring non-conformal cells.", 1); - res = Staircaser{ res }.getSelectiveMesh(nonConformalCells,Staircaser::GapsFillingType::Split); + res = Staircaser{ res }.getSelectiveMesh(nonConformalCells,Staircaser::GapsFillingType::Insert); RedundancyCleaner::removeOverlappedDimensionOneAndLowerElementsAndEquivalentSurfaces(res); logNumberOfTriangles(countMeshElementsIf(res, isTriangle)); diff --git a/tmp/sphere.tessellator.json b/tmp/sphere.tessellator.json new file mode 100644 index 0000000..8601a88 --- /dev/null +++ b/tmp/sphere.tessellator.json @@ -0,0 +1,10 @@ +{ + "grid": { + "numberOfCells": [50, 50, 50], + "boundingBox": [ + [-100.0, -100.0, -100.0], + [ 100.0, 100.0, 100.0] + ] + }, + "object": {"filename": "sphere.stl"} +} \ No newline at end of file From 000c861525a7aa8a82f613fb6b86b7faa1537973 Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Fri, 31 Oct 2025 11:54:53 +0100 Subject: [PATCH 06/23] Minor update --- testData/cases/sphere/sphere.conformal.tessellator.json | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testData/cases/sphere/sphere.conformal.tessellator.json b/testData/cases/sphere/sphere.conformal.tessellator.json index 67d7ef2..c32cc53 100644 --- a/testData/cases/sphere/sphere.conformal.tessellator.json +++ b/testData/cases/sphere/sphere.conformal.tessellator.json @@ -1,6 +1,6 @@ { "grid": { - "numberOfCells": [50, 50, 50], + "numberOfCells": [30, 30, 30], "boundingBox": [ [-100.0, -100.0, -100.0], [ 100.0, 100.0, 100.0] From b694c5b9aeb019cdc830bbf7b04db59227999dbb Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Mon, 3 Nov 2025 13:49:38 +0100 Subject: [PATCH 07/23] Refactor findCommonNeighborsVertices to handle missing cell elements and improve boundary coordinate indexing in getSelectiveMesh --- src/core/Staircaser.cpp | 60 ++++++++++++++++++++++++++--------------- 1 file changed, 38 insertions(+), 22 deletions(-) diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index 7ee6ff4..f921dd1 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -72,31 +72,41 @@ IdSet findCommonNeighborsVertices(const Mesh& mesh, const std::pairvertices; - auto it = std::find(elementVertices.begin(), elementVertices.end(), vertex1); - - if (it != elementVertices.end()) { - auto pos = std::distance(elementVertices.begin(), it); - auto n = elementVertices.size(); - - neighborsOfVertex1.insert(elementVertices[(pos + n - 1) % n]); - neighborsOfVertex1.insert(elementVertices[(pos + 1) % n]); - } + auto it = cellElemMap.find(c); + if (it == cellElemMap.end()) { + continue; + } else { + for (const auto& e : cellElemMap.at(c)) { + const auto& elementVertices = e->vertices; + auto it = std::find(elementVertices.begin(), elementVertices.end(), vertex1); + + if (it != elementVertices.end()) { + auto pos = std::distance(elementVertices.begin(), it); + auto n = elementVertices.size(); + + neighborsOfVertex1.insert(elementVertices[(pos + n - 1) % n]); + neighborsOfVertex1.insert(elementVertices[(pos + 1) % n]); + } + } } } for (const auto& c : cellsCoord2) { - for (const auto& e : cellElemMap.at(c)) { - const auto& elementVertices = e->vertices; - auto it = std::find(elementVertices.begin(), elementVertices.end(), vertex2); - - if (it != elementVertices.end()) { - auto pos = std::distance(elementVertices.begin(), it); - auto n = elementVertices.size(); - - neighborsOfVertex2.insert(elementVertices[(pos + n - 1) % n]); - neighborsOfVertex2.insert(elementVertices[(pos + 1) % n]); + auto it = cellElemMap.find(c); + if (it == cellElemMap.end()) { + continue; + } else { + for (const auto& e : cellElemMap.at(c)) { + const auto& elementVertices = e->vertices; + auto it = std::find(elementVertices.begin(), elementVertices.end(), vertex2); + + if (it != elementVertices.end()) { + auto pos = std::distance(elementVertices.begin(), it); + auto n = elementVertices.size(); + + neighborsOfVertex2.insert(elementVertices[(pos + n - 1) % n]); + neighborsOfVertex2.insert(elementVertices[(pos + 1) % n]); + } } } } @@ -185,7 +195,13 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi CoordinateId newIndex; if (isOnCellBoundary) { auto it = coordinateMap.find(toRelative(calculateStaircasedCell(vertexCoord))); - newIndex = it->second; + if (it == coordinateMap.end()) { + mesh_.coordinates.push_back(toRelative(calculateStaircasedCell(vertexCoord))); + newIndex = int (mesh_.coordinates.size() - 1); + coordinateMap.emplace(toRelative(calculateStaircasedCell(vertexCoord)), newIndex); + } else { + newIndex = it->second; + } boundaryCoordinates.push_back(toRelative(calculateStaircasedCell(vertexCoord))); } else { auto it = coordinateMap.find(vertexCoord); From fc39c70b9995b7cf7bbe81ecb41c81335d2fe048 Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Tue, 4 Nov 2025 12:37:52 +0100 Subject: [PATCH 08/23] Refactor getSelectiveMesh to improve boundary checking for structured cells --- src/core/Staircaser.cpp | 36 ++++++++++++++++++++++++------------ 1 file changed, 24 insertions(+), 12 deletions(-) diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index f921dd1..0d4c3b2 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -216,27 +216,39 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi newElement.vertices.push_back(newIndex); } - bool isAllCoordinatesOnCellBoundary = true; + bool isAllCoordinatesOnTheSameCellBoundary = false; + std::set commonStructuredCells; + bool firstVertex = true; for (const auto& vertex : newElement.vertices) { const auto& vertexCoord = mesh_.coordinates[vertex]; auto touchingCells = GridTools::getTouchingCells(vertexCoord); - - bool vertexOnBoundary = false; - for (const auto& touchingCell : touchingCells) { - if (cellsToStructure.count(touchingCell)) { - vertexOnBoundary = true; - break; - } + + std::set structuredTouchingCells; + for (const auto& c : touchingCells) { + if (cellsToStructure.count(c)) { + structuredTouchingCells.insert(c); + } } - if(!vertexOnBoundary) { - isAllCoordinatesOnCellBoundary = false; - break; + if (firstVertex) { + commonStructuredCells = std::move(structuredTouchingCells); + firstVertex = false; + } else { + std::set intersection; + std::set_intersection(commonStructuredCells.begin(), commonStructuredCells.end(), + structuredTouchingCells.begin(), structuredTouchingCells.end(), + std::inserter(intersection, intersection.begin())); + commonStructuredCells = std::move(intersection); } + + if (commonStructuredCells.empty()) break; } - if (!isAllCoordinatesOnCellBoundary) { + isAllCoordinatesOnTheSameCellBoundary = !commonStructuredCells.empty(); + + + if (!isAllCoordinatesOnTheSameCellBoundary) { meshGroup.elements.push_back(newElement); for (size_t i = 0; i < boundaryCoordinates.size(); ++i) { From d7d79bf1575b0212bbd2e10253fb00e911701748 Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Thu, 6 Nov 2025 14:26:47 +0100 Subject: [PATCH 09/23] Enhance getSelectiveMesh to support spliting triangles when a surface is structured to a line --- src/core/Staircaser.cpp | 125 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 125 insertions(+) diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index 0d4c3b2..363d774 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -168,6 +168,8 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi CoordinateMap coordinateMap = buildCoordinateMap(mesh_.coordinates); meshGroup.elements.reserve(meshGroup.elements.size() + inputGroup.elements.size()); + std::map> cellElemMap_withNewElements; + for (const auto& [cell, elements] : cellElemMap) { if (cellsToStructure.count(cell) ) { continue; @@ -250,6 +252,7 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi if (!isAllCoordinatesOnTheSameCellBoundary) { meshGroup.elements.push_back(newElement); + cellElemMap_withNewElements[cell].push_back(&meshGroup.elements.back()); for (size_t i = 0; i < boundaryCoordinates.size(); ++i) { for (size_t j = i + 1; j < boundaryCoordinates.size(); ++j) { @@ -258,6 +261,128 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi } } } + + // Now lets split some elements according to the starcased cells + std::set elementsConvertedInLines; + auto itCell = cellElemMap_withNewElements.find(cell); + if (itCell != cellElemMap_withNewElements.end()) { + + for(const auto& e: cellElemMap_withNewElements.at(cell)) { + if (!e->vertices.empty()) { + bool allVerticesAreEqual = std::all_of( + e->vertices.begin() + 1, + e->vertices.end(), + [&](int v) { return v == e->vertices.front(); } + ); + + if (allVerticesAreEqual) { + continue; + } + } + + if (e->isTriangle()) { + const auto& firstVertexCoords = mesh_.coordinates[e->vertices[0]]; + const auto& secondVertexCoords = mesh_.coordinates[e->vertices[1]]; + const auto& thirdVertexCoords = mesh_.coordinates[e->vertices[2]]; + + int equalCoords = 0; + + bool sameAxis = false; + for (int dim = 0; dim < 3; ++dim) { + if (firstVertexCoords[dim] == secondVertexCoords[dim] && + secondVertexCoords[dim] == thirdVertexCoords[dim]) { + ++equalCoords; + } + } + + bool isTwoVerticesEqual = (firstVertexCoords == secondVertexCoords) || + (secondVertexCoords == thirdVertexCoords) || + (firstVertexCoords == thirdVertexCoords); + + if(equalCoords == 2 && !isTwoVerticesEqual) { + elementsConvertedInLines.insert(*e); + } + } + } + + for (const auto& e: elementsConvertedInLines) { + for (const auto& otherElementsInCell: cellElemMap_withNewElements.at(cell)) { + if (e == *otherElementsInCell) { + continue; + } + + IdSet sharedVertices; + int sharedCount = 0; + int lastCommonVertexIndex = -2; + bool correctOrientation = false; + + for (int i = 0; i < otherElementsInCell->vertices.size(); ++i) { + int vOther = otherElementsInCell->vertices[i]; + + if (std::find(e.vertices.begin(), e.vertices.end(), vOther) != e.vertices.end()) { + sharedVertices.insert(vOther); + ++sharedCount; + } + + if (i == lastCommonVertexIndex + 1) { + correctOrientation = true; + } + } + + bool shareExactlyTwo = (sharedCount == 2); + + if(shareExactlyTwo && otherElementsInCell->isTriangle()) { + + auto v1 = *sharedVertices.begin(); + auto v2 = *std::next(sharedVertices.begin()); + + if(!correctOrientation) { + std::swap(v1, v2); + } + + auto itV1 = std::find(e.vertices.begin(), e.vertices.end(), v1); + auto itV2 = std::find(otherElementsInCell->vertices.begin(), otherElementsInCell->vertices.end(), v2); + auto itV4 = e.vertices.end(); + auto itV3 = otherElementsInCell->vertices.end(); + + if(itV1 == e.vertices.begin()) { + itV4 = e.vertices.end() - 1; + } else { + itV4 = std::prev(itV1); + } + if(itV2 == otherElementsInCell->vertices.begin()) { + itV3 = otherElementsInCell->vertices.end() - 1; + } else { + itV3 = std::prev(itV2); + } + + const auto& v3 = *itV3; + const auto& v4 = *itV4; + + Element triangle1; + triangle1.type = Element::Type::Surface; + triangle1.vertices = { v3, v2, v4 }; + + Element triangle2; + triangle2.type = Element::Type::Surface; + triangle2.vertices = { v4, v1, v3 }; + + meshGroup.elements.push_back(triangle1); + meshGroup.elements.push_back(triangle2); + + auto itOtherElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), *otherElementsInCell); + auto itElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), e); + + if (itOtherElement != meshGroup.elements.end()) { + meshGroup.elements.erase(itOtherElement); + } + if (itElement != meshGroup.elements.end()) { + meshGroup.elements.erase(itElement); + } + } + } + } + } } } From 98a7b38e9da91162e128bad369d9306af0728371 Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Thu, 6 Nov 2025 15:35:29 +0100 Subject: [PATCH 10/23] Refactor getSelectiveMesh to improve vertex comparison logic and streamline triangle processing --- src/core/Staircaser.cpp | 37 +++++++++++++++++++------------------ 1 file changed, 19 insertions(+), 18 deletions(-) diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index 363d774..c69bbe4 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -295,17 +295,18 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi } } - bool isTwoVerticesEqual = (firstVertexCoords == secondVertexCoords) || + bool areTwoVerticesEqual = (firstVertexCoords == secondVertexCoords) || (secondVertexCoords == thirdVertexCoords) || (firstVertexCoords == thirdVertexCoords); - if(equalCoords == 2 && !isTwoVerticesEqual) { + if(equalCoords == 2 && !areTwoVerticesEqual) { elementsConvertedInLines.insert(*e); } } } for (const auto& e: elementsConvertedInLines) { + bool processed = false; for (const auto& otherElementsInCell: cellElemMap_withNewElements.at(cell)) { if (e == *otherElementsInCell) { continue; @@ -340,25 +341,22 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi std::swap(v1, v2); } - auto itV1 = std::find(e.vertices.begin(), e.vertices.end(), v1); - auto itV2 = std::find(otherElementsInCell->vertices.begin(), otherElementsInCell->vertices.end(), v2); - auto itV4 = e.vertices.end(); - auto itV3 = otherElementsInCell->vertices.end(); - - if(itV1 == e.vertices.begin()) { - itV4 = e.vertices.end() - 1; - } else { - itV4 = std::prev(itV1); + CoordinateId v4 = -1; + for (const auto& v : e.vertices) { + if (v != v1 && v != v2) { + v4 = v; + break; + } } - if(itV2 == otherElementsInCell->vertices.begin()) { - itV3 = otherElementsInCell->vertices.end() - 1; - } else { - itV3 = std::prev(itV2); + + CoordinateId v3 = -1; + for (const auto& v : otherElementsInCell->vertices) { + if (v != v1 && v != v2) { + v3 = v; + break; + } } - const auto& v3 = *itV3; - const auto& v4 = *itV4; - Element triangle1; triangle1.type = Element::Type::Surface; triangle1.vertices = { v3, v2, v4 }; @@ -379,8 +377,11 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi if (itElement != meshGroup.elements.end()) { meshGroup.elements.erase(itElement); } + processed = true; + break; } } + if(processed) {continue;} } } } From 84f901d1df8c5858b5cc1a4f25ffe89a342e2919 Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Fri, 7 Nov 2025 14:55:00 +0100 Subject: [PATCH 11/23] Refactor getSelectiveMesh to optimize element removal and enhance redundancy cleaning --- src/core/Staircaser.cpp | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index c69bbe4..b9c9707 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -143,6 +143,7 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi fillerType_ = type; RelativePairSet boundaryCoordinatePairs; + std::vector> toRemove(mesh_.groups.size()); for (std::size_t g = 0; g < mesh_.groups.size(); ++g) { auto& inputGroup = inputMesh_.groups[g]; @@ -364,19 +365,25 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi Element triangle2; triangle2.type = Element::Type::Surface; triangle2.vertices = { v4, v1, v3 }; - - meshGroup.elements.push_back(triangle1); - meshGroup.elements.push_back(triangle2); - + auto itOtherElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), *otherElementsInCell); - auto itElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), e); - if (itOtherElement != meshGroup.elements.end()) { - meshGroup.elements.erase(itOtherElement); + ElementId otherElementId = std::distance(meshGroup.elements.begin(), itOtherElement); + toRemove[g].insert(otherElementId); } + + auto itElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), e); if (itElement != meshGroup.elements.end()) { - meshGroup.elements.erase(itElement); + ElementId elementId = std::distance(meshGroup.elements.begin(), itElement); + toRemove[g].insert(elementId); } + + meshGroup.elements.push_back(triangle1); + meshGroup.elements.push_back(triangle2); + + RedundancyCleaner::removeElements(mesh_, toRemove); + toRemove[g].clear(); + processed = true; break; } @@ -386,9 +393,10 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi } } } - + RedundancyCleaner::fuseCoords(mesh_); RedundancyCleaner::removeDegenerateElements(mesh_); + RedundancyCleaner::removeRepeatedElements(mesh_); RedundancyCleaner::cleanCoords(mesh_); for (auto it = boundaryCoordinatePairs.begin(); it != boundaryCoordinatePairs.end();) { From 6791dbe29358d42acd672b2af7537b4303fc72b9 Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Fri, 7 Nov 2025 14:55:21 +0100 Subject: [PATCH 12/23] Add test for selective structurer to split lines with neighboring triangles into two another triangles --- test/core/StaircaserTest.cpp | 102 +++++++++++++++++++++++++++++++++++ 1 file changed, 102 insertions(+) diff --git a/test/core/StaircaserTest.cpp b/test/core/StaircaserTest.cpp index 8820542..6c8540b 100644 --- a/test/core/StaircaserTest.cpp +++ b/test/core/StaircaserTest.cpp @@ -2641,3 +2641,105 @@ TEST_F(StaircaserTest, selectiveStructurerFillingGapsInFrontier_Insert) } } } + +TEST_F(StaircaserTest, selectiveStructurer_SplitLinesWithNeighborTriangle) +{ + // *----0========2=============5 *---(0->3)==(2->5)========(5->6) + // | ‾-_ ║|\ ║ | ‾-_ ║ \ ║ + // | ‾-_ ║ | \ ║ | ‾-_ ║ \ ║ + // | - ║ | \ ║ | - ║ \ ║ + // | ‾1 | \ ║ | (1->4) \ ║ + // | |‾- | \ ║ | ║‾‾--__ \ ║ + // | | ‾-_| \ ║ | ║ ‾‾--_\ ║ + // *-------------*-----3=======4 -> *-----------(3->0)========(4->1) + // | | \ | | | ║ + // | | \ | | | ║ + // | | \ | | | ║ + // | | \ | | | ║ + // | | \ | | | ║ + // | | \ | | | ║ + // *-------------*-------------6 *-------------*-----------(6->2) + + float lowerCoordinateValue = -5.0; + float upperCoordinateValue = 5.0; + int numberOfCells = 3; + float step = 5.0; + assert((upperCoordinateValue - lowerCoordinateValue) / (numberOfCells - 1) == step); + + std::set cellSet; + cellSet.insert(Cell({1, 0, 0})); + + Mesh mesh; + mesh.grid = GridTools::buildCartesianGrid(lowerCoordinateValue, upperCoordinateValue, numberOfCells); + mesh.coordinates = { + Relative({ 0.4, 2.0, 1.0 }), // 0 + Relative({ 1.0, 1.4, 1.0 }), // 1 + Relative({ 1.0, 2.0, 1.0 }), // 2 + Relative({ 1.4, 1.0, 1.0 }), // 3 + Relative({ 2.0, 1.0, 1.0 }), // 4 + Relative({ 2.0, 2.0, 1.0 }), // 5 + Relative({ 2.0, 0.0, 1.0 }), // 6 + }; + + mesh.groups.resize(1); + mesh.groups[0].elements = { + Element({0, 1, 2}, Element::Type::Surface), + Element({1, 3, 2}, Element::Type::Surface), + Element({3, 4, 2}, Element::Type::Surface), + Element({2, 4, 5}, Element::Type::Surface), + Element({3, 6, 4}, Element::Type::Surface), + }; + + Relatives expectedRelatives = { + Relative({ 1.0, 1.0, 1.0 }), // (3->0) + Relative({ 2.0, 1.0, 1.0 }), // (4->1) + Relative({ 2.0, 0.0, 1.0 }), // (6->2) + Relative({ 0.4, 2.0, 1.0 }), // (0->3) + Relative({ 1.0, 1.4, 1.0 }), // (1->4) + Relative({ 1.0, 2.0, 1.0 }), // (2->5) + Relative({ 2.0, 2.0, 1.0 }), // (5->6) + }; + + Elements expectedElements = { + Element({0, 1}, Element::Type::Line), + Element({1, 2}, Element::Type::Line), + Element({2, 1}, Element::Type::Line), + Element({1, 0}, Element::Type::Line), + Element({3, 4, 5}, Element::Type::Surface), + Element({5, 1, 6}, Element::Type::Surface), + Element({1, 0, 4}, Element::Type::Surface), + Element({4, 5, 1}, Element::Type::Surface), + }; + + auto resultMesh = Staircaser{ mesh }.getSelectiveMesh(cellSet); + + ASSERT_EQ(resultMesh.coordinates.size(), expectedRelatives.size()); + ASSERT_EQ(resultMesh.groups.size(), 1); + ASSERT_EQ(resultMesh.groups[0].elements.size(), expectedElements.size()); + + for (std::size_t i = 0; i < resultMesh.coordinates.size(); ++i) { + for (std::size_t axis = 0; axis < 3; ++axis) { + EXPECT_EQ(resultMesh.coordinates[i][axis], expectedRelatives[i][axis]); + } + } + + ASSERT_TRUE(resultMesh.groups[0].elements[0].isLine()); + ASSERT_TRUE(resultMesh.groups[0].elements[1].isLine()); + ASSERT_TRUE(resultMesh.groups[0].elements[2].isLine()); + ASSERT_TRUE(resultMesh.groups[0].elements[3].isLine()); + + ASSERT_TRUE(resultMesh.groups[0].elements[4].isTriangle()); + ASSERT_TRUE(resultMesh.groups[0].elements[5].isTriangle()); + ASSERT_TRUE(resultMesh.groups[0].elements[6].isTriangle()); + ASSERT_TRUE(resultMesh.groups[0].elements[7].isTriangle()); + + for (std::size_t e = 0; e < expectedElements.size(); ++e) { + auto& resultElement = resultMesh.groups[0].elements[e]; + auto& expectedElement = expectedElements[e]; + + for (std::size_t v = 0; v < expectedElement.vertices.size(); ++v) { + EXPECT_EQ(resultElement.vertices[v], expectedElement.vertices[v]); + } + } + +} From 0124dbcbccea5e547a1bd2e5af72ebbf87682047 Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Mon, 10 Nov 2025 09:21:20 +0100 Subject: [PATCH 13/23] Refactor getSelectiveMesh to correct vertex orientation logic and update test for selective structurer to reflect changes in triangle definitions --- src/core/Staircaser.cpp | 10 +++++----- test/core/StaircaserTest.cpp | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index b9c9707..5b4737d 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -315,7 +315,6 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi IdSet sharedVertices; int sharedCount = 0; - int lastCommonVertexIndex = -2; bool correctOrientation = false; for (int i = 0; i < otherElementsInCell->vertices.size(); ++i) { @@ -325,10 +324,7 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi sharedVertices.insert(vOther); ++sharedCount; } - - if (i == lastCommonVertexIndex + 1) { - correctOrientation = true; - } + } bool shareExactlyTwo = (sharedCount == 2); @@ -337,6 +333,10 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi auto v1 = *sharedVertices.begin(); auto v2 = *std::next(sharedVertices.begin()); + + auto itV1 = std::find(e.vertices.begin(), e.vertices.end(), v1); + auto itV2 = std::find(e.vertices.begin(), e.vertices.end(), v2); + correctOrientation = (itV1 < itV2); if(!correctOrientation) { std::swap(v1, v2); diff --git a/test/core/StaircaserTest.cpp b/test/core/StaircaserTest.cpp index 6c8540b..7191d27 100644 --- a/test/core/StaircaserTest.cpp +++ b/test/core/StaircaserTest.cpp @@ -2707,8 +2707,8 @@ TEST_F(StaircaserTest, selectiveStructurer_SplitLinesWithNeighborTriangle) Element({1, 0}, Element::Type::Line), Element({3, 4, 5}, Element::Type::Surface), Element({5, 1, 6}, Element::Type::Surface), - Element({1, 0, 4}, Element::Type::Surface), - Element({4, 5, 1}, Element::Type::Surface), + Element({1, 5, 4}, Element::Type::Surface), + Element({4, 0, 1}, Element::Type::Surface), }; auto resultMesh = Staircaser{ mesh }.getSelectiveMesh(cellSet); From 27096404457b03e12c3e3f6be99cbc825c15b296 Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Mon, 10 Nov 2025 12:17:33 +0100 Subject: [PATCH 14/23] Refactor getSelectiveMesh to streamline element processing and introduce helper functions for better readability and maintainability --- src/core/Staircaser.cpp | 421 +++++++++++++++++++++------------------- src/core/Staircaser.h | 4 + 2 files changed, 225 insertions(+), 200 deletions(-) diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index 5b4737d..5d09b17 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -166,92 +166,20 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi } } - CoordinateMap coordinateMap = buildCoordinateMap(mesh_.coordinates); - meshGroup.elements.reserve(meshGroup.elements.size() + inputGroup.elements.size()); std::map> cellElemMap_withNewElements; + CoordinateMap coordinateMap = buildCoordinateMap(mesh_.coordinates); for (const auto& [cell, elements] : cellElemMap) { if (cellsToStructure.count(cell) ) { continue; } for (const auto e : elements) { - Element newElement; - newElement.type = e->type; - newElement.vertices.reserve(e->vertices.size()); - - Relatives boundaryCoordinates; - - for (const auto& vertexIndex : e->vertices) { - const auto& vertexCoord = inputMesh_.coordinates[vertexIndex]; - - bool isOnCellBoundary = false; - auto touchingCells = GridTools::getTouchingCells(vertexCoord); - - for (const auto& touchingCell : touchingCells) { - if (cellsToStructure.count(touchingCell)) { - isOnCellBoundary = true; - break; - } - } - - CoordinateId newIndex; - if (isOnCellBoundary) { - auto it = coordinateMap.find(toRelative(calculateStaircasedCell(vertexCoord))); - if (it == coordinateMap.end()) { - mesh_.coordinates.push_back(toRelative(calculateStaircasedCell(vertexCoord))); - newIndex = int (mesh_.coordinates.size() - 1); - coordinateMap.emplace(toRelative(calculateStaircasedCell(vertexCoord)), newIndex); - } else { - newIndex = it->second; - } - boundaryCoordinates.push_back(toRelative(calculateStaircasedCell(vertexCoord))); - } else { - auto it = coordinateMap.find(vertexCoord); - if (it == coordinateMap.end()) { - mesh_.coordinates.push_back(vertexCoord); - newIndex = int (mesh_.coordinates.size() - 1); - coordinateMap.emplace(vertexCoord, newIndex); - } else { - newIndex = it->second; - } - } - newElement.vertices.push_back(newIndex); - } - - bool isAllCoordinatesOnTheSameCellBoundary = false; - std::set commonStructuredCells; - bool firstVertex = true; - - for (const auto& vertex : newElement.vertices) { - const auto& vertexCoord = mesh_.coordinates[vertex]; - auto touchingCells = GridTools::getTouchingCells(vertexCoord); - - std::set structuredTouchingCells; - for (const auto& c : touchingCells) { - if (cellsToStructure.count(c)) { - structuredTouchingCells.insert(c); - } - } - - if (firstVertex) { - commonStructuredCells = std::move(structuredTouchingCells); - firstVertex = false; - } else { - std::set intersection; - std::set_intersection(commonStructuredCells.begin(), commonStructuredCells.end(), - structuredTouchingCells.begin(), structuredTouchingCells.end(), - std::inserter(intersection, intersection.begin())); - commonStructuredCells = std::move(intersection); - } - - if (commonStructuredCells.empty()) break; - } - - isAllCoordinatesOnTheSameCellBoundary = !commonStructuredCells.empty(); + auto [newElement, boundaryCoordinates] = obtainNewIndexForElement(*e, cellsToStructure, coordinateMap); + auto allCoordsOnBoundary = isAllCoordinatesOnTheSameCellBoundary(newElement, cellsToStructure); - if (!isAllCoordinatesOnTheSameCellBoundary) { + if (!allCoordsOnBoundary) { meshGroup.elements.push_back(newElement); cellElemMap_withNewElements[cell].push_back(&meshGroup.elements.back()); @@ -263,133 +191,12 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi } } - // Now lets split some elements according to the starcased cells - std::set elementsConvertedInLines; auto itCell = cellElemMap_withNewElements.find(cell); if (itCell != cellElemMap_withNewElements.end()) { - for(const auto& e: cellElemMap_withNewElements.at(cell)) { - if (!e->vertices.empty()) { - bool allVerticesAreEqual = std::all_of( - e->vertices.begin() + 1, - e->vertices.end(), - [&](int v) { return v == e->vertices.front(); } - ); - - if (allVerticesAreEqual) { - continue; - } - } - - if (e->isTriangle()) { - const auto& firstVertexCoords = mesh_.coordinates[e->vertices[0]]; - const auto& secondVertexCoords = mesh_.coordinates[e->vertices[1]]; - const auto& thirdVertexCoords = mesh_.coordinates[e->vertices[2]]; - - int equalCoords = 0; - - bool sameAxis = false; - for (int dim = 0; dim < 3; ++dim) { - if (firstVertexCoords[dim] == secondVertexCoords[dim] && - secondVertexCoords[dim] == thirdVertexCoords[dim]) { - ++equalCoords; - } - } - - bool areTwoVerticesEqual = (firstVertexCoords == secondVertexCoords) || - (secondVertexCoords == thirdVertexCoords) || - (firstVertexCoords == thirdVertexCoords); - - if(equalCoords == 2 && !areTwoVerticesEqual) { - elementsConvertedInLines.insert(*e); - } - } - } - - for (const auto& e: elementsConvertedInLines) { - bool processed = false; - for (const auto& otherElementsInCell: cellElemMap_withNewElements.at(cell)) { - if (e == *otherElementsInCell) { - continue; - } - - IdSet sharedVertices; - int sharedCount = 0; - bool correctOrientation = false; - - for (int i = 0; i < otherElementsInCell->vertices.size(); ++i) { - int vOther = otherElementsInCell->vertices[i]; - - if (std::find(e.vertices.begin(), e.vertices.end(), vOther) != e.vertices.end()) { - sharedVertices.insert(vOther); - ++sharedCount; - } - - } - - bool shareExactlyTwo = (sharedCount == 2); - - if(shareExactlyTwo && otherElementsInCell->isTriangle()) { - - auto v1 = *sharedVertices.begin(); - auto v2 = *std::next(sharedVertices.begin()); - - auto itV1 = std::find(e.vertices.begin(), e.vertices.end(), v1); - auto itV2 = std::find(e.vertices.begin(), e.vertices.end(), v2); - correctOrientation = (itV1 < itV2); - - if(!correctOrientation) { - std::swap(v1, v2); - } - - CoordinateId v4 = -1; - for (const auto& v : e.vertices) { - if (v != v1 && v != v2) { - v4 = v; - break; - } - } - - CoordinateId v3 = -1; - for (const auto& v : otherElementsInCell->vertices) { - if (v != v1 && v != v2) { - v3 = v; - break; - } - } - - Element triangle1; - triangle1.type = Element::Type::Surface; - triangle1.vertices = { v3, v2, v4 }; - - Element triangle2; - triangle2.type = Element::Type::Surface; - triangle2.vertices = { v4, v1, v3 }; - - auto itOtherElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), *otherElementsInCell); - if (itOtherElement != meshGroup.elements.end()) { - ElementId otherElementId = std::distance(meshGroup.elements.begin(), itOtherElement); - toRemove[g].insert(otherElementId); - } - - auto itElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), e); - if (itElement != meshGroup.elements.end()) { - ElementId elementId = std::distance(meshGroup.elements.begin(), itElement); - toRemove[g].insert(elementId); - } - - meshGroup.elements.push_back(triangle1); - meshGroup.elements.push_back(triangle2); - - RedundancyCleaner::removeElements(mesh_, toRemove); - toRemove[g].clear(); - - processed = true; - break; - } - } - if(processed) {continue;} - } + auto elementsConvertedInLines = getElementsConvertedInLines(cell, cellElemMap_withNewElements); + splitLinesWithNeighborTriangle(g, elementsConvertedInLines, meshGroup, cell, cellElemMap_withNewElements, toRemove); + } } } @@ -555,6 +362,220 @@ void Staircaser::fillGaps(const RelativePairSet boundaryCoordinatePairs) { RedundancyCleaner::removeDegenerateElements(mesh_); } +std::pair Staircaser::obtainNewIndexForElement(const Element& e, const std::set& cellsToStructure, CoordinateMap& coordinateMap) { + Element newElement; + newElement.type = e.type; + newElement.vertices.reserve(e.vertices.size()); + + Relatives boundaryCoordinates; + + for (const auto& vertexIndex : e.vertices) { + const auto& vertexCoord = inputMesh_.coordinates[vertexIndex]; + + bool isOnCellBoundary = false; + auto touchingCells = GridTools::getTouchingCells(vertexCoord); + + for (const auto& touchingCell : touchingCells) { + if (cellsToStructure.count(touchingCell)) { + isOnCellBoundary = true; + break; + } + } + + CoordinateId newIndex; + if (isOnCellBoundary) { + auto relCoord = toRelative(calculateStaircasedCell(vertexCoord)); + auto it = coordinateMap.find(relCoord); + + if (it == coordinateMap.end()) { + mesh_.coordinates.push_back(relCoord); + newIndex = int(mesh_.coordinates.size() - 1); + coordinateMap.emplace(relCoord, newIndex); + } else { + newIndex = it->second; + } + + boundaryCoordinates.push_back(relCoord); + } else { + auto it = coordinateMap.find(vertexCoord); + if (it == coordinateMap.end()) { + mesh_.coordinates.push_back(vertexCoord); + newIndex = int(mesh_.coordinates.size() - 1); + coordinateMap.emplace(vertexCoord, newIndex); + } else { + newIndex = it->second; + } + } + + newElement.vertices.push_back(newIndex); + } + + return {newElement, boundaryCoordinates}; +} + +bool Staircaser::isAllCoordinatesOnTheSameCellBoundary(const Element& newElement, const std::set& cellsToStructure) { + std::set commonStructuredCells; + bool firstVertex = true; + + for (const auto& vertex : newElement.vertices) { + const auto& vertexCoord = mesh_.coordinates[vertex]; + auto touchingCells = GridTools::getTouchingCells(vertexCoord); + + std::set structuredTouchingCells; + for (const auto& c : touchingCells) { + if (cellsToStructure.count(c)) { + structuredTouchingCells.insert(c); + } + } + + if (firstVertex) { + commonStructuredCells = std::move(structuredTouchingCells); + firstVertex = false; + } else { + std::set intersection; + std::set_intersection(commonStructuredCells.begin(), commonStructuredCells.end(), + structuredTouchingCells.begin(), structuredTouchingCells.end(), + std::inserter(intersection, intersection.begin())); + commonStructuredCells = std::move(intersection); + } + + if (commonStructuredCells.empty()) break; + } + + return !commonStructuredCells.empty(); +} + +std::set Staircaser::getElementsConvertedInLines(Cell cell, std::map> cellElemMap) { + std::set elementsConvertedInLines; + + for (const auto& e : cellElemMap.at(cell)) { + if (!e->vertices.empty()) { + bool allVerticesAreEqual = std::all_of( + e->vertices.begin() + 1, + e->vertices.end(), + [&](int v) { return v == e->vertices.front(); } + ); + + if (allVerticesAreEqual) { + continue; + } + } + + if (e->isTriangle()) { + const auto& firstVertexCoords = mesh_.coordinates[e->vertices[0]]; + const auto& secondVertexCoords = mesh_.coordinates[e->vertices[1]]; + const auto& thirdVertexCoords = mesh_.coordinates[e->vertices[2]]; + + int equalCoords = 0; + + for (int dim = 0; dim < 3; ++dim) { + if (firstVertexCoords[dim] == secondVertexCoords[dim] && + secondVertexCoords[dim] == thirdVertexCoords[dim]) { + ++equalCoords; + } + } + + bool areTwoVerticesEqual = + (firstVertexCoords == secondVertexCoords) || + (secondVertexCoords == thirdVertexCoords) || + (firstVertexCoords == thirdVertexCoords); + + if (equalCoords == 2 && !areTwoVerticesEqual) { + elementsConvertedInLines.insert(*e); + } + } + } + + return elementsConvertedInLines; +} + +void Staircaser::splitLinesWithNeighborTriangle(const size_t groupIndex, std::set elementsConvertedInLines, Group& meshGroup, const Cell cell, std::map> cellElemMap, std::vector> toRemove) { + for (const auto& e: elementsConvertedInLines) { + bool processed = false; + for (const auto& otherElementsInCell: cellElemMap.at(cell)) { + if (e == *otherElementsInCell) { + continue; + } + + IdSet sharedVertices; + int sharedCount = 0; + bool correctOrientation = false; + + for (int i = 0; i < otherElementsInCell->vertices.size(); ++i) { + int vOther = otherElementsInCell->vertices[i]; + + if (std::find(e.vertices.begin(), e.vertices.end(), vOther) != e.vertices.end()) { + sharedVertices.insert(vOther); + ++sharedCount; + } + + } + + bool shareExactlyTwo = (sharedCount == 2); + + if(shareExactlyTwo && otherElementsInCell->isTriangle()) { + + auto v1 = *sharedVertices.begin(); + auto v2 = *std::next(sharedVertices.begin()); + + auto itV1 = std::find(e.vertices.begin(), e.vertices.end(), v1); + auto itV2 = std::find(e.vertices.begin(), e.vertices.end(), v2); + correctOrientation = (itV1 < itV2); + + if(!correctOrientation) { + std::swap(v1, v2); + } + + CoordinateId v4 = -1; + for (const auto& v : e.vertices) { + if (v != v1 && v != v2) { + v4 = v; + break; + } + } + + CoordinateId v3 = -1; + for (const auto& v : otherElementsInCell->vertices) { + if (v != v1 && v != v2) { + v3 = v; + break; + } + } + + Element triangle1; + triangle1.type = Element::Type::Surface; + triangle1.vertices = { v3, v2, v4 }; + + Element triangle2; + triangle2.type = Element::Type::Surface; + triangle2.vertices = { v4, v1, v3 }; + + auto itOtherElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), *otherElementsInCell); + if (itOtherElement != meshGroup.elements.end()) { + ElementId otherElementId = std::distance(meshGroup.elements.begin(), itOtherElement); + toRemove[groupIndex].insert(otherElementId); + } + + auto itElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), e); + if (itElement != meshGroup.elements.end()) { + ElementId elementId = std::distance(meshGroup.elements.begin(), itElement); + toRemove[groupIndex].insert(elementId); + } + + meshGroup.elements.push_back(triangle1); + meshGroup.elements.push_back(triangle2); + + RedundancyCleaner::removeElements(mesh_, toRemove); + toRemove[groupIndex].clear(); + + processed = true; + break; + } + } + if(processed) {continue;} + } +} + void Staircaser::processTriangleAndAddToGroup(const Element& triangle, const Relatives& originalRelatives, Group& group){ Group edges; diff --git a/src/core/Staircaser.h b/src/core/Staircaser.h index 77220e9..b4a4d20 100644 --- a/src/core/Staircaser.h +++ b/src/core/Staircaser.h @@ -70,6 +70,10 @@ class Staircaser : public utils::GridTools { std::vector calculateMiddleCellsBetweenTwoRelatives(Relative& startExtreme, Relative& endExtreme); void calculateRelativeIdSetByCellSurface(const Relatives& relatives, std::map& relativesByCellSurface); void fillGaps(const RelativePairSet boundaryCoordinatePairs); + std::pair obtainNewIndexForElement(const Element& e, const std::set& cellsToStructure, CoordinateMap& coordinateMap); + bool isAllCoordinatesOnTheSameCellBoundary(const Element& newElement, const std::set& cellsToStructure); + std::set getElementsConvertedInLines(const Cell cell, std::map> cellElemMap); + void splitLinesWithNeighborTriangle(const size_t groupIndex, std::set elementsConvertedInLines, Group& meshGroup, const Cell cell, std::map> cellElemMap, std::vector> toRemove); }; From 2cfb09522c0d5c9054a3464acd7d4d889e096a43 Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Mon, 3 Nov 2025 13:49:38 +0100 Subject: [PATCH 15/23] Refactor findCommonNeighborsVertices to handle missing cell elements and improve boundary coordinate indexing in getSelectiveMesh --- src/core/Staircaser.cpp | 60 ++++++++++++++++++++++++++--------------- 1 file changed, 38 insertions(+), 22 deletions(-) diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index 7ee6ff4..f921dd1 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -72,31 +72,41 @@ IdSet findCommonNeighborsVertices(const Mesh& mesh, const std::pairvertices; - auto it = std::find(elementVertices.begin(), elementVertices.end(), vertex1); - - if (it != elementVertices.end()) { - auto pos = std::distance(elementVertices.begin(), it); - auto n = elementVertices.size(); - - neighborsOfVertex1.insert(elementVertices[(pos + n - 1) % n]); - neighborsOfVertex1.insert(elementVertices[(pos + 1) % n]); - } + auto it = cellElemMap.find(c); + if (it == cellElemMap.end()) { + continue; + } else { + for (const auto& e : cellElemMap.at(c)) { + const auto& elementVertices = e->vertices; + auto it = std::find(elementVertices.begin(), elementVertices.end(), vertex1); + + if (it != elementVertices.end()) { + auto pos = std::distance(elementVertices.begin(), it); + auto n = elementVertices.size(); + + neighborsOfVertex1.insert(elementVertices[(pos + n - 1) % n]); + neighborsOfVertex1.insert(elementVertices[(pos + 1) % n]); + } + } } } for (const auto& c : cellsCoord2) { - for (const auto& e : cellElemMap.at(c)) { - const auto& elementVertices = e->vertices; - auto it = std::find(elementVertices.begin(), elementVertices.end(), vertex2); - - if (it != elementVertices.end()) { - auto pos = std::distance(elementVertices.begin(), it); - auto n = elementVertices.size(); - - neighborsOfVertex2.insert(elementVertices[(pos + n - 1) % n]); - neighborsOfVertex2.insert(elementVertices[(pos + 1) % n]); + auto it = cellElemMap.find(c); + if (it == cellElemMap.end()) { + continue; + } else { + for (const auto& e : cellElemMap.at(c)) { + const auto& elementVertices = e->vertices; + auto it = std::find(elementVertices.begin(), elementVertices.end(), vertex2); + + if (it != elementVertices.end()) { + auto pos = std::distance(elementVertices.begin(), it); + auto n = elementVertices.size(); + + neighborsOfVertex2.insert(elementVertices[(pos + n - 1) % n]); + neighborsOfVertex2.insert(elementVertices[(pos + 1) % n]); + } } } } @@ -185,7 +195,13 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi CoordinateId newIndex; if (isOnCellBoundary) { auto it = coordinateMap.find(toRelative(calculateStaircasedCell(vertexCoord))); - newIndex = it->second; + if (it == coordinateMap.end()) { + mesh_.coordinates.push_back(toRelative(calculateStaircasedCell(vertexCoord))); + newIndex = int (mesh_.coordinates.size() - 1); + coordinateMap.emplace(toRelative(calculateStaircasedCell(vertexCoord)), newIndex); + } else { + newIndex = it->second; + } boundaryCoordinates.push_back(toRelative(calculateStaircasedCell(vertexCoord))); } else { auto it = coordinateMap.find(vertexCoord); From 6427a795782dc51f15d14481da589f4cca7d8b89 Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Tue, 4 Nov 2025 12:37:52 +0100 Subject: [PATCH 16/23] Refactor getSelectiveMesh to improve boundary checking for structured cells --- src/core/Staircaser.cpp | 36 ++++++++++++++++++++++++------------ 1 file changed, 24 insertions(+), 12 deletions(-) diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index f921dd1..0d4c3b2 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -216,27 +216,39 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi newElement.vertices.push_back(newIndex); } - bool isAllCoordinatesOnCellBoundary = true; + bool isAllCoordinatesOnTheSameCellBoundary = false; + std::set commonStructuredCells; + bool firstVertex = true; for (const auto& vertex : newElement.vertices) { const auto& vertexCoord = mesh_.coordinates[vertex]; auto touchingCells = GridTools::getTouchingCells(vertexCoord); - - bool vertexOnBoundary = false; - for (const auto& touchingCell : touchingCells) { - if (cellsToStructure.count(touchingCell)) { - vertexOnBoundary = true; - break; - } + + std::set structuredTouchingCells; + for (const auto& c : touchingCells) { + if (cellsToStructure.count(c)) { + structuredTouchingCells.insert(c); + } } - if(!vertexOnBoundary) { - isAllCoordinatesOnCellBoundary = false; - break; + if (firstVertex) { + commonStructuredCells = std::move(structuredTouchingCells); + firstVertex = false; + } else { + std::set intersection; + std::set_intersection(commonStructuredCells.begin(), commonStructuredCells.end(), + structuredTouchingCells.begin(), structuredTouchingCells.end(), + std::inserter(intersection, intersection.begin())); + commonStructuredCells = std::move(intersection); } + + if (commonStructuredCells.empty()) break; } - if (!isAllCoordinatesOnCellBoundary) { + isAllCoordinatesOnTheSameCellBoundary = !commonStructuredCells.empty(); + + + if (!isAllCoordinatesOnTheSameCellBoundary) { meshGroup.elements.push_back(newElement); for (size_t i = 0; i < boundaryCoordinates.size(); ++i) { From b713cbdfd6dca0b8dd14d30e4e6e0b5c3e128616 Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Thu, 6 Nov 2025 14:26:47 +0100 Subject: [PATCH 17/23] Enhance getSelectiveMesh to support spliting triangles when a surface is structured to a line --- src/core/Staircaser.cpp | 125 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 125 insertions(+) diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index 0d4c3b2..363d774 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -168,6 +168,8 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi CoordinateMap coordinateMap = buildCoordinateMap(mesh_.coordinates); meshGroup.elements.reserve(meshGroup.elements.size() + inputGroup.elements.size()); + std::map> cellElemMap_withNewElements; + for (const auto& [cell, elements] : cellElemMap) { if (cellsToStructure.count(cell) ) { continue; @@ -250,6 +252,7 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi if (!isAllCoordinatesOnTheSameCellBoundary) { meshGroup.elements.push_back(newElement); + cellElemMap_withNewElements[cell].push_back(&meshGroup.elements.back()); for (size_t i = 0; i < boundaryCoordinates.size(); ++i) { for (size_t j = i + 1; j < boundaryCoordinates.size(); ++j) { @@ -258,6 +261,128 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi } } } + + // Now lets split some elements according to the starcased cells + std::set elementsConvertedInLines; + auto itCell = cellElemMap_withNewElements.find(cell); + if (itCell != cellElemMap_withNewElements.end()) { + + for(const auto& e: cellElemMap_withNewElements.at(cell)) { + if (!e->vertices.empty()) { + bool allVerticesAreEqual = std::all_of( + e->vertices.begin() + 1, + e->vertices.end(), + [&](int v) { return v == e->vertices.front(); } + ); + + if (allVerticesAreEqual) { + continue; + } + } + + if (e->isTriangle()) { + const auto& firstVertexCoords = mesh_.coordinates[e->vertices[0]]; + const auto& secondVertexCoords = mesh_.coordinates[e->vertices[1]]; + const auto& thirdVertexCoords = mesh_.coordinates[e->vertices[2]]; + + int equalCoords = 0; + + bool sameAxis = false; + for (int dim = 0; dim < 3; ++dim) { + if (firstVertexCoords[dim] == secondVertexCoords[dim] && + secondVertexCoords[dim] == thirdVertexCoords[dim]) { + ++equalCoords; + } + } + + bool isTwoVerticesEqual = (firstVertexCoords == secondVertexCoords) || + (secondVertexCoords == thirdVertexCoords) || + (firstVertexCoords == thirdVertexCoords); + + if(equalCoords == 2 && !isTwoVerticesEqual) { + elementsConvertedInLines.insert(*e); + } + } + } + + for (const auto& e: elementsConvertedInLines) { + for (const auto& otherElementsInCell: cellElemMap_withNewElements.at(cell)) { + if (e == *otherElementsInCell) { + continue; + } + + IdSet sharedVertices; + int sharedCount = 0; + int lastCommonVertexIndex = -2; + bool correctOrientation = false; + + for (int i = 0; i < otherElementsInCell->vertices.size(); ++i) { + int vOther = otherElementsInCell->vertices[i]; + + if (std::find(e.vertices.begin(), e.vertices.end(), vOther) != e.vertices.end()) { + sharedVertices.insert(vOther); + ++sharedCount; + } + + if (i == lastCommonVertexIndex + 1) { + correctOrientation = true; + } + } + + bool shareExactlyTwo = (sharedCount == 2); + + if(shareExactlyTwo && otherElementsInCell->isTriangle()) { + + auto v1 = *sharedVertices.begin(); + auto v2 = *std::next(sharedVertices.begin()); + + if(!correctOrientation) { + std::swap(v1, v2); + } + + auto itV1 = std::find(e.vertices.begin(), e.vertices.end(), v1); + auto itV2 = std::find(otherElementsInCell->vertices.begin(), otherElementsInCell->vertices.end(), v2); + auto itV4 = e.vertices.end(); + auto itV3 = otherElementsInCell->vertices.end(); + + if(itV1 == e.vertices.begin()) { + itV4 = e.vertices.end() - 1; + } else { + itV4 = std::prev(itV1); + } + if(itV2 == otherElementsInCell->vertices.begin()) { + itV3 = otherElementsInCell->vertices.end() - 1; + } else { + itV3 = std::prev(itV2); + } + + const auto& v3 = *itV3; + const auto& v4 = *itV4; + + Element triangle1; + triangle1.type = Element::Type::Surface; + triangle1.vertices = { v3, v2, v4 }; + + Element triangle2; + triangle2.type = Element::Type::Surface; + triangle2.vertices = { v4, v1, v3 }; + + meshGroup.elements.push_back(triangle1); + meshGroup.elements.push_back(triangle2); + + auto itOtherElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), *otherElementsInCell); + auto itElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), e); + + if (itOtherElement != meshGroup.elements.end()) { + meshGroup.elements.erase(itOtherElement); + } + if (itElement != meshGroup.elements.end()) { + meshGroup.elements.erase(itElement); + } + } + } + } + } } } From 9738333da39503d3ea5af5f148092ddfc1a60301 Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Thu, 6 Nov 2025 15:35:29 +0100 Subject: [PATCH 18/23] Refactor getSelectiveMesh to improve vertex comparison logic and streamline triangle processing --- src/core/Staircaser.cpp | 37 +++++++++++++++++++------------------ 1 file changed, 19 insertions(+), 18 deletions(-) diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index 363d774..c69bbe4 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -295,17 +295,18 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi } } - bool isTwoVerticesEqual = (firstVertexCoords == secondVertexCoords) || + bool areTwoVerticesEqual = (firstVertexCoords == secondVertexCoords) || (secondVertexCoords == thirdVertexCoords) || (firstVertexCoords == thirdVertexCoords); - if(equalCoords == 2 && !isTwoVerticesEqual) { + if(equalCoords == 2 && !areTwoVerticesEqual) { elementsConvertedInLines.insert(*e); } } } for (const auto& e: elementsConvertedInLines) { + bool processed = false; for (const auto& otherElementsInCell: cellElemMap_withNewElements.at(cell)) { if (e == *otherElementsInCell) { continue; @@ -340,25 +341,22 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi std::swap(v1, v2); } - auto itV1 = std::find(e.vertices.begin(), e.vertices.end(), v1); - auto itV2 = std::find(otherElementsInCell->vertices.begin(), otherElementsInCell->vertices.end(), v2); - auto itV4 = e.vertices.end(); - auto itV3 = otherElementsInCell->vertices.end(); - - if(itV1 == e.vertices.begin()) { - itV4 = e.vertices.end() - 1; - } else { - itV4 = std::prev(itV1); + CoordinateId v4 = -1; + for (const auto& v : e.vertices) { + if (v != v1 && v != v2) { + v4 = v; + break; + } } - if(itV2 == otherElementsInCell->vertices.begin()) { - itV3 = otherElementsInCell->vertices.end() - 1; - } else { - itV3 = std::prev(itV2); + + CoordinateId v3 = -1; + for (const auto& v : otherElementsInCell->vertices) { + if (v != v1 && v != v2) { + v3 = v; + break; + } } - const auto& v3 = *itV3; - const auto& v4 = *itV4; - Element triangle1; triangle1.type = Element::Type::Surface; triangle1.vertices = { v3, v2, v4 }; @@ -379,8 +377,11 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi if (itElement != meshGroup.elements.end()) { meshGroup.elements.erase(itElement); } + processed = true; + break; } } + if(processed) {continue;} } } } From 5233c8f9673b8b7fcb2cdead2893fa1eee3508bc Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Fri, 7 Nov 2025 14:55:00 +0100 Subject: [PATCH 19/23] Refactor getSelectiveMesh to optimize element removal and enhance redundancy cleaning --- src/core/Staircaser.cpp | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index c69bbe4..b9c9707 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -143,6 +143,7 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi fillerType_ = type; RelativePairSet boundaryCoordinatePairs; + std::vector> toRemove(mesh_.groups.size()); for (std::size_t g = 0; g < mesh_.groups.size(); ++g) { auto& inputGroup = inputMesh_.groups[g]; @@ -364,19 +365,25 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi Element triangle2; triangle2.type = Element::Type::Surface; triangle2.vertices = { v4, v1, v3 }; - - meshGroup.elements.push_back(triangle1); - meshGroup.elements.push_back(triangle2); - + auto itOtherElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), *otherElementsInCell); - auto itElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), e); - if (itOtherElement != meshGroup.elements.end()) { - meshGroup.elements.erase(itOtherElement); + ElementId otherElementId = std::distance(meshGroup.elements.begin(), itOtherElement); + toRemove[g].insert(otherElementId); } + + auto itElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), e); if (itElement != meshGroup.elements.end()) { - meshGroup.elements.erase(itElement); + ElementId elementId = std::distance(meshGroup.elements.begin(), itElement); + toRemove[g].insert(elementId); } + + meshGroup.elements.push_back(triangle1); + meshGroup.elements.push_back(triangle2); + + RedundancyCleaner::removeElements(mesh_, toRemove); + toRemove[g].clear(); + processed = true; break; } @@ -386,9 +393,10 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi } } } - + RedundancyCleaner::fuseCoords(mesh_); RedundancyCleaner::removeDegenerateElements(mesh_); + RedundancyCleaner::removeRepeatedElements(mesh_); RedundancyCleaner::cleanCoords(mesh_); for (auto it = boundaryCoordinatePairs.begin(); it != boundaryCoordinatePairs.end();) { From ee66b39dc1db22f3f8418a8e859417b95b690f7a Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Fri, 7 Nov 2025 14:55:21 +0100 Subject: [PATCH 20/23] Add test for selective structurer to split lines with neighboring triangles into two another triangles --- test/core/StaircaserTest.cpp | 102 +++++++++++++++++++++++++++++++++++ 1 file changed, 102 insertions(+) diff --git a/test/core/StaircaserTest.cpp b/test/core/StaircaserTest.cpp index 8820542..6c8540b 100644 --- a/test/core/StaircaserTest.cpp +++ b/test/core/StaircaserTest.cpp @@ -2641,3 +2641,105 @@ TEST_F(StaircaserTest, selectiveStructurerFillingGapsInFrontier_Insert) } } } + +TEST_F(StaircaserTest, selectiveStructurer_SplitLinesWithNeighborTriangle) +{ + // *----0========2=============5 *---(0->3)==(2->5)========(5->6) + // | ‾-_ ║|\ ║ | ‾-_ ║ \ ║ + // | ‾-_ ║ | \ ║ | ‾-_ ║ \ ║ + // | - ║ | \ ║ | - ║ \ ║ + // | ‾1 | \ ║ | (1->4) \ ║ + // | |‾- | \ ║ | ║‾‾--__ \ ║ + // | | ‾-_| \ ║ | ║ ‾‾--_\ ║ + // *-------------*-----3=======4 -> *-----------(3->0)========(4->1) + // | | \ | | | ║ + // | | \ | | | ║ + // | | \ | | | ║ + // | | \ | | | ║ + // | | \ | | | ║ + // | | \ | | | ║ + // *-------------*-------------6 *-------------*-----------(6->2) + + float lowerCoordinateValue = -5.0; + float upperCoordinateValue = 5.0; + int numberOfCells = 3; + float step = 5.0; + assert((upperCoordinateValue - lowerCoordinateValue) / (numberOfCells - 1) == step); + + std::set cellSet; + cellSet.insert(Cell({1, 0, 0})); + + Mesh mesh; + mesh.grid = GridTools::buildCartesianGrid(lowerCoordinateValue, upperCoordinateValue, numberOfCells); + mesh.coordinates = { + Relative({ 0.4, 2.0, 1.0 }), // 0 + Relative({ 1.0, 1.4, 1.0 }), // 1 + Relative({ 1.0, 2.0, 1.0 }), // 2 + Relative({ 1.4, 1.0, 1.0 }), // 3 + Relative({ 2.0, 1.0, 1.0 }), // 4 + Relative({ 2.0, 2.0, 1.0 }), // 5 + Relative({ 2.0, 0.0, 1.0 }), // 6 + }; + + mesh.groups.resize(1); + mesh.groups[0].elements = { + Element({0, 1, 2}, Element::Type::Surface), + Element({1, 3, 2}, Element::Type::Surface), + Element({3, 4, 2}, Element::Type::Surface), + Element({2, 4, 5}, Element::Type::Surface), + Element({3, 6, 4}, Element::Type::Surface), + }; + + Relatives expectedRelatives = { + Relative({ 1.0, 1.0, 1.0 }), // (3->0) + Relative({ 2.0, 1.0, 1.0 }), // (4->1) + Relative({ 2.0, 0.0, 1.0 }), // (6->2) + Relative({ 0.4, 2.0, 1.0 }), // (0->3) + Relative({ 1.0, 1.4, 1.0 }), // (1->4) + Relative({ 1.0, 2.0, 1.0 }), // (2->5) + Relative({ 2.0, 2.0, 1.0 }), // (5->6) + }; + + Elements expectedElements = { + Element({0, 1}, Element::Type::Line), + Element({1, 2}, Element::Type::Line), + Element({2, 1}, Element::Type::Line), + Element({1, 0}, Element::Type::Line), + Element({3, 4, 5}, Element::Type::Surface), + Element({5, 1, 6}, Element::Type::Surface), + Element({1, 0, 4}, Element::Type::Surface), + Element({4, 5, 1}, Element::Type::Surface), + }; + + auto resultMesh = Staircaser{ mesh }.getSelectiveMesh(cellSet); + + ASSERT_EQ(resultMesh.coordinates.size(), expectedRelatives.size()); + ASSERT_EQ(resultMesh.groups.size(), 1); + ASSERT_EQ(resultMesh.groups[0].elements.size(), expectedElements.size()); + + for (std::size_t i = 0; i < resultMesh.coordinates.size(); ++i) { + for (std::size_t axis = 0; axis < 3; ++axis) { + EXPECT_EQ(resultMesh.coordinates[i][axis], expectedRelatives[i][axis]); + } + } + + ASSERT_TRUE(resultMesh.groups[0].elements[0].isLine()); + ASSERT_TRUE(resultMesh.groups[0].elements[1].isLine()); + ASSERT_TRUE(resultMesh.groups[0].elements[2].isLine()); + ASSERT_TRUE(resultMesh.groups[0].elements[3].isLine()); + + ASSERT_TRUE(resultMesh.groups[0].elements[4].isTriangle()); + ASSERT_TRUE(resultMesh.groups[0].elements[5].isTriangle()); + ASSERT_TRUE(resultMesh.groups[0].elements[6].isTriangle()); + ASSERT_TRUE(resultMesh.groups[0].elements[7].isTriangle()); + + for (std::size_t e = 0; e < expectedElements.size(); ++e) { + auto& resultElement = resultMesh.groups[0].elements[e]; + auto& expectedElement = expectedElements[e]; + + for (std::size_t v = 0; v < expectedElement.vertices.size(); ++v) { + EXPECT_EQ(resultElement.vertices[v], expectedElement.vertices[v]); + } + } + +} From 93c43732de023771d4b2a66f7d0503adc5a1eea5 Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Mon, 10 Nov 2025 09:21:20 +0100 Subject: [PATCH 21/23] Refactor getSelectiveMesh to correct vertex orientation logic and update test for selective structurer to reflect changes in triangle definitions --- src/core/Staircaser.cpp | 10 +++++----- test/core/StaircaserTest.cpp | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index b9c9707..5b4737d 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -315,7 +315,6 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi IdSet sharedVertices; int sharedCount = 0; - int lastCommonVertexIndex = -2; bool correctOrientation = false; for (int i = 0; i < otherElementsInCell->vertices.size(); ++i) { @@ -325,10 +324,7 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi sharedVertices.insert(vOther); ++sharedCount; } - - if (i == lastCommonVertexIndex + 1) { - correctOrientation = true; - } + } bool shareExactlyTwo = (sharedCount == 2); @@ -337,6 +333,10 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi auto v1 = *sharedVertices.begin(); auto v2 = *std::next(sharedVertices.begin()); + + auto itV1 = std::find(e.vertices.begin(), e.vertices.end(), v1); + auto itV2 = std::find(e.vertices.begin(), e.vertices.end(), v2); + correctOrientation = (itV1 < itV2); if(!correctOrientation) { std::swap(v1, v2); diff --git a/test/core/StaircaserTest.cpp b/test/core/StaircaserTest.cpp index 6c8540b..7191d27 100644 --- a/test/core/StaircaserTest.cpp +++ b/test/core/StaircaserTest.cpp @@ -2707,8 +2707,8 @@ TEST_F(StaircaserTest, selectiveStructurer_SplitLinesWithNeighborTriangle) Element({1, 0}, Element::Type::Line), Element({3, 4, 5}, Element::Type::Surface), Element({5, 1, 6}, Element::Type::Surface), - Element({1, 0, 4}, Element::Type::Surface), - Element({4, 5, 1}, Element::Type::Surface), + Element({1, 5, 4}, Element::Type::Surface), + Element({4, 0, 1}, Element::Type::Surface), }; auto resultMesh = Staircaser{ mesh }.getSelectiveMesh(cellSet); From cbc941aabe82f866516ba4371cc88b58145fea85 Mon Sep 17 00:00:00 2001 From: ashybabashyba Date: Mon, 10 Nov 2025 12:17:33 +0100 Subject: [PATCH 22/23] Refactor getSelectiveMesh to streamline element processing and introduce helper functions for better readability and maintainability --- src/core/Staircaser.cpp | 421 +++++++++++++++++++++------------------- src/core/Staircaser.h | 4 + 2 files changed, 225 insertions(+), 200 deletions(-) diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index 5b4737d..5d09b17 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -166,92 +166,20 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi } } - CoordinateMap coordinateMap = buildCoordinateMap(mesh_.coordinates); - meshGroup.elements.reserve(meshGroup.elements.size() + inputGroup.elements.size()); std::map> cellElemMap_withNewElements; + CoordinateMap coordinateMap = buildCoordinateMap(mesh_.coordinates); for (const auto& [cell, elements] : cellElemMap) { if (cellsToStructure.count(cell) ) { continue; } for (const auto e : elements) { - Element newElement; - newElement.type = e->type; - newElement.vertices.reserve(e->vertices.size()); - - Relatives boundaryCoordinates; - - for (const auto& vertexIndex : e->vertices) { - const auto& vertexCoord = inputMesh_.coordinates[vertexIndex]; - - bool isOnCellBoundary = false; - auto touchingCells = GridTools::getTouchingCells(vertexCoord); - - for (const auto& touchingCell : touchingCells) { - if (cellsToStructure.count(touchingCell)) { - isOnCellBoundary = true; - break; - } - } - - CoordinateId newIndex; - if (isOnCellBoundary) { - auto it = coordinateMap.find(toRelative(calculateStaircasedCell(vertexCoord))); - if (it == coordinateMap.end()) { - mesh_.coordinates.push_back(toRelative(calculateStaircasedCell(vertexCoord))); - newIndex = int (mesh_.coordinates.size() - 1); - coordinateMap.emplace(toRelative(calculateStaircasedCell(vertexCoord)), newIndex); - } else { - newIndex = it->second; - } - boundaryCoordinates.push_back(toRelative(calculateStaircasedCell(vertexCoord))); - } else { - auto it = coordinateMap.find(vertexCoord); - if (it == coordinateMap.end()) { - mesh_.coordinates.push_back(vertexCoord); - newIndex = int (mesh_.coordinates.size() - 1); - coordinateMap.emplace(vertexCoord, newIndex); - } else { - newIndex = it->second; - } - } - newElement.vertices.push_back(newIndex); - } - - bool isAllCoordinatesOnTheSameCellBoundary = false; - std::set commonStructuredCells; - bool firstVertex = true; - - for (const auto& vertex : newElement.vertices) { - const auto& vertexCoord = mesh_.coordinates[vertex]; - auto touchingCells = GridTools::getTouchingCells(vertexCoord); - - std::set structuredTouchingCells; - for (const auto& c : touchingCells) { - if (cellsToStructure.count(c)) { - structuredTouchingCells.insert(c); - } - } - - if (firstVertex) { - commonStructuredCells = std::move(structuredTouchingCells); - firstVertex = false; - } else { - std::set intersection; - std::set_intersection(commonStructuredCells.begin(), commonStructuredCells.end(), - structuredTouchingCells.begin(), structuredTouchingCells.end(), - std::inserter(intersection, intersection.begin())); - commonStructuredCells = std::move(intersection); - } - - if (commonStructuredCells.empty()) break; - } - - isAllCoordinatesOnTheSameCellBoundary = !commonStructuredCells.empty(); + auto [newElement, boundaryCoordinates] = obtainNewIndexForElement(*e, cellsToStructure, coordinateMap); + auto allCoordsOnBoundary = isAllCoordinatesOnTheSameCellBoundary(newElement, cellsToStructure); - if (!isAllCoordinatesOnTheSameCellBoundary) { + if (!allCoordsOnBoundary) { meshGroup.elements.push_back(newElement); cellElemMap_withNewElements[cell].push_back(&meshGroup.elements.back()); @@ -263,133 +191,12 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi } } - // Now lets split some elements according to the starcased cells - std::set elementsConvertedInLines; auto itCell = cellElemMap_withNewElements.find(cell); if (itCell != cellElemMap_withNewElements.end()) { - for(const auto& e: cellElemMap_withNewElements.at(cell)) { - if (!e->vertices.empty()) { - bool allVerticesAreEqual = std::all_of( - e->vertices.begin() + 1, - e->vertices.end(), - [&](int v) { return v == e->vertices.front(); } - ); - - if (allVerticesAreEqual) { - continue; - } - } - - if (e->isTriangle()) { - const auto& firstVertexCoords = mesh_.coordinates[e->vertices[0]]; - const auto& secondVertexCoords = mesh_.coordinates[e->vertices[1]]; - const auto& thirdVertexCoords = mesh_.coordinates[e->vertices[2]]; - - int equalCoords = 0; - - bool sameAxis = false; - for (int dim = 0; dim < 3; ++dim) { - if (firstVertexCoords[dim] == secondVertexCoords[dim] && - secondVertexCoords[dim] == thirdVertexCoords[dim]) { - ++equalCoords; - } - } - - bool areTwoVerticesEqual = (firstVertexCoords == secondVertexCoords) || - (secondVertexCoords == thirdVertexCoords) || - (firstVertexCoords == thirdVertexCoords); - - if(equalCoords == 2 && !areTwoVerticesEqual) { - elementsConvertedInLines.insert(*e); - } - } - } - - for (const auto& e: elementsConvertedInLines) { - bool processed = false; - for (const auto& otherElementsInCell: cellElemMap_withNewElements.at(cell)) { - if (e == *otherElementsInCell) { - continue; - } - - IdSet sharedVertices; - int sharedCount = 0; - bool correctOrientation = false; - - for (int i = 0; i < otherElementsInCell->vertices.size(); ++i) { - int vOther = otherElementsInCell->vertices[i]; - - if (std::find(e.vertices.begin(), e.vertices.end(), vOther) != e.vertices.end()) { - sharedVertices.insert(vOther); - ++sharedCount; - } - - } - - bool shareExactlyTwo = (sharedCount == 2); - - if(shareExactlyTwo && otherElementsInCell->isTriangle()) { - - auto v1 = *sharedVertices.begin(); - auto v2 = *std::next(sharedVertices.begin()); - - auto itV1 = std::find(e.vertices.begin(), e.vertices.end(), v1); - auto itV2 = std::find(e.vertices.begin(), e.vertices.end(), v2); - correctOrientation = (itV1 < itV2); - - if(!correctOrientation) { - std::swap(v1, v2); - } - - CoordinateId v4 = -1; - for (const auto& v : e.vertices) { - if (v != v1 && v != v2) { - v4 = v; - break; - } - } - - CoordinateId v3 = -1; - for (const auto& v : otherElementsInCell->vertices) { - if (v != v1 && v != v2) { - v3 = v; - break; - } - } - - Element triangle1; - triangle1.type = Element::Type::Surface; - triangle1.vertices = { v3, v2, v4 }; - - Element triangle2; - triangle2.type = Element::Type::Surface; - triangle2.vertices = { v4, v1, v3 }; - - auto itOtherElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), *otherElementsInCell); - if (itOtherElement != meshGroup.elements.end()) { - ElementId otherElementId = std::distance(meshGroup.elements.begin(), itOtherElement); - toRemove[g].insert(otherElementId); - } - - auto itElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), e); - if (itElement != meshGroup.elements.end()) { - ElementId elementId = std::distance(meshGroup.elements.begin(), itElement); - toRemove[g].insert(elementId); - } - - meshGroup.elements.push_back(triangle1); - meshGroup.elements.push_back(triangle2); - - RedundancyCleaner::removeElements(mesh_, toRemove); - toRemove[g].clear(); - - processed = true; - break; - } - } - if(processed) {continue;} - } + auto elementsConvertedInLines = getElementsConvertedInLines(cell, cellElemMap_withNewElements); + splitLinesWithNeighborTriangle(g, elementsConvertedInLines, meshGroup, cell, cellElemMap_withNewElements, toRemove); + } } } @@ -555,6 +362,220 @@ void Staircaser::fillGaps(const RelativePairSet boundaryCoordinatePairs) { RedundancyCleaner::removeDegenerateElements(mesh_); } +std::pair Staircaser::obtainNewIndexForElement(const Element& e, const std::set& cellsToStructure, CoordinateMap& coordinateMap) { + Element newElement; + newElement.type = e.type; + newElement.vertices.reserve(e.vertices.size()); + + Relatives boundaryCoordinates; + + for (const auto& vertexIndex : e.vertices) { + const auto& vertexCoord = inputMesh_.coordinates[vertexIndex]; + + bool isOnCellBoundary = false; + auto touchingCells = GridTools::getTouchingCells(vertexCoord); + + for (const auto& touchingCell : touchingCells) { + if (cellsToStructure.count(touchingCell)) { + isOnCellBoundary = true; + break; + } + } + + CoordinateId newIndex; + if (isOnCellBoundary) { + auto relCoord = toRelative(calculateStaircasedCell(vertexCoord)); + auto it = coordinateMap.find(relCoord); + + if (it == coordinateMap.end()) { + mesh_.coordinates.push_back(relCoord); + newIndex = int(mesh_.coordinates.size() - 1); + coordinateMap.emplace(relCoord, newIndex); + } else { + newIndex = it->second; + } + + boundaryCoordinates.push_back(relCoord); + } else { + auto it = coordinateMap.find(vertexCoord); + if (it == coordinateMap.end()) { + mesh_.coordinates.push_back(vertexCoord); + newIndex = int(mesh_.coordinates.size() - 1); + coordinateMap.emplace(vertexCoord, newIndex); + } else { + newIndex = it->second; + } + } + + newElement.vertices.push_back(newIndex); + } + + return {newElement, boundaryCoordinates}; +} + +bool Staircaser::isAllCoordinatesOnTheSameCellBoundary(const Element& newElement, const std::set& cellsToStructure) { + std::set commonStructuredCells; + bool firstVertex = true; + + for (const auto& vertex : newElement.vertices) { + const auto& vertexCoord = mesh_.coordinates[vertex]; + auto touchingCells = GridTools::getTouchingCells(vertexCoord); + + std::set structuredTouchingCells; + for (const auto& c : touchingCells) { + if (cellsToStructure.count(c)) { + structuredTouchingCells.insert(c); + } + } + + if (firstVertex) { + commonStructuredCells = std::move(structuredTouchingCells); + firstVertex = false; + } else { + std::set intersection; + std::set_intersection(commonStructuredCells.begin(), commonStructuredCells.end(), + structuredTouchingCells.begin(), structuredTouchingCells.end(), + std::inserter(intersection, intersection.begin())); + commonStructuredCells = std::move(intersection); + } + + if (commonStructuredCells.empty()) break; + } + + return !commonStructuredCells.empty(); +} + +std::set Staircaser::getElementsConvertedInLines(Cell cell, std::map> cellElemMap) { + std::set elementsConvertedInLines; + + for (const auto& e : cellElemMap.at(cell)) { + if (!e->vertices.empty()) { + bool allVerticesAreEqual = std::all_of( + e->vertices.begin() + 1, + e->vertices.end(), + [&](int v) { return v == e->vertices.front(); } + ); + + if (allVerticesAreEqual) { + continue; + } + } + + if (e->isTriangle()) { + const auto& firstVertexCoords = mesh_.coordinates[e->vertices[0]]; + const auto& secondVertexCoords = mesh_.coordinates[e->vertices[1]]; + const auto& thirdVertexCoords = mesh_.coordinates[e->vertices[2]]; + + int equalCoords = 0; + + for (int dim = 0; dim < 3; ++dim) { + if (firstVertexCoords[dim] == secondVertexCoords[dim] && + secondVertexCoords[dim] == thirdVertexCoords[dim]) { + ++equalCoords; + } + } + + bool areTwoVerticesEqual = + (firstVertexCoords == secondVertexCoords) || + (secondVertexCoords == thirdVertexCoords) || + (firstVertexCoords == thirdVertexCoords); + + if (equalCoords == 2 && !areTwoVerticesEqual) { + elementsConvertedInLines.insert(*e); + } + } + } + + return elementsConvertedInLines; +} + +void Staircaser::splitLinesWithNeighborTriangle(const size_t groupIndex, std::set elementsConvertedInLines, Group& meshGroup, const Cell cell, std::map> cellElemMap, std::vector> toRemove) { + for (const auto& e: elementsConvertedInLines) { + bool processed = false; + for (const auto& otherElementsInCell: cellElemMap.at(cell)) { + if (e == *otherElementsInCell) { + continue; + } + + IdSet sharedVertices; + int sharedCount = 0; + bool correctOrientation = false; + + for (int i = 0; i < otherElementsInCell->vertices.size(); ++i) { + int vOther = otherElementsInCell->vertices[i]; + + if (std::find(e.vertices.begin(), e.vertices.end(), vOther) != e.vertices.end()) { + sharedVertices.insert(vOther); + ++sharedCount; + } + + } + + bool shareExactlyTwo = (sharedCount == 2); + + if(shareExactlyTwo && otherElementsInCell->isTriangle()) { + + auto v1 = *sharedVertices.begin(); + auto v2 = *std::next(sharedVertices.begin()); + + auto itV1 = std::find(e.vertices.begin(), e.vertices.end(), v1); + auto itV2 = std::find(e.vertices.begin(), e.vertices.end(), v2); + correctOrientation = (itV1 < itV2); + + if(!correctOrientation) { + std::swap(v1, v2); + } + + CoordinateId v4 = -1; + for (const auto& v : e.vertices) { + if (v != v1 && v != v2) { + v4 = v; + break; + } + } + + CoordinateId v3 = -1; + for (const auto& v : otherElementsInCell->vertices) { + if (v != v1 && v != v2) { + v3 = v; + break; + } + } + + Element triangle1; + triangle1.type = Element::Type::Surface; + triangle1.vertices = { v3, v2, v4 }; + + Element triangle2; + triangle2.type = Element::Type::Surface; + triangle2.vertices = { v4, v1, v3 }; + + auto itOtherElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), *otherElementsInCell); + if (itOtherElement != meshGroup.elements.end()) { + ElementId otherElementId = std::distance(meshGroup.elements.begin(), itOtherElement); + toRemove[groupIndex].insert(otherElementId); + } + + auto itElement = std::find(meshGroup.elements.begin(), meshGroup.elements.end(), e); + if (itElement != meshGroup.elements.end()) { + ElementId elementId = std::distance(meshGroup.elements.begin(), itElement); + toRemove[groupIndex].insert(elementId); + } + + meshGroup.elements.push_back(triangle1); + meshGroup.elements.push_back(triangle2); + + RedundancyCleaner::removeElements(mesh_, toRemove); + toRemove[groupIndex].clear(); + + processed = true; + break; + } + } + if(processed) {continue;} + } +} + void Staircaser::processTriangleAndAddToGroup(const Element& triangle, const Relatives& originalRelatives, Group& group){ Group edges; diff --git a/src/core/Staircaser.h b/src/core/Staircaser.h index 77220e9..b4a4d20 100644 --- a/src/core/Staircaser.h +++ b/src/core/Staircaser.h @@ -70,6 +70,10 @@ class Staircaser : public utils::GridTools { std::vector calculateMiddleCellsBetweenTwoRelatives(Relative& startExtreme, Relative& endExtreme); void calculateRelativeIdSetByCellSurface(const Relatives& relatives, std::map& relativesByCellSurface); void fillGaps(const RelativePairSet boundaryCoordinatePairs); + std::pair obtainNewIndexForElement(const Element& e, const std::set& cellsToStructure, CoordinateMap& coordinateMap); + bool isAllCoordinatesOnTheSameCellBoundary(const Element& newElement, const std::set& cellsToStructure); + std::set getElementsConvertedInLines(const Cell cell, std::map> cellElemMap); + void splitLinesWithNeighborTriangle(const size_t groupIndex, std::set elementsConvertedInLines, Group& meshGroup, const Cell cell, std::map> cellElemMap, std::vector> toRemove); }; From 5ac9e7aea74687b6fa43cc7fae76bee0b6922247 Mon Sep 17 00:00:00 2001 From: Alberto-o Date: Tue, 11 Nov 2025 12:48:18 +0100 Subject: [PATCH 23/23] Minor --- src/app/launcher.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/app/launcher.h b/src/app/launcher.h index ca5388b..e54b2f5 100644 --- a/src/app/launcher.h +++ b/src/app/launcher.h @@ -6,7 +6,7 @@ namespace meshlib::app { const std::string conformal_mesher ("conformal"); -const std::string structured_mesher ("structured"); +const std::string staircase_mesher ("staircase"); int launcher(int argc, const char* argv[]);