diff --git a/.github/workflows/build-and-test.yml b/.github/workflows/build-and-test.yml index c6b8acd..94cfdac 100644 --- a/.github/workflows/build-and-test.yml +++ b/.github/workflows/build-and-test.yml @@ -28,8 +28,8 @@ jobs: name: ${{ matrix.preset.os }} / ${{matrix.preset.name}} / ${{matrix.build-type}} runs-on: ${{ matrix.preset.os }} - env: - VCPKG_BINARY_SOURCES: clear + # Inherit the workflow-level x-gha binary cache. Do not override this + # with `clear`: Windows VTK is otherwise rebuilt from source on every run. steps: - name: checkout repository diff --git a/.gitignore b/.gitignore index 130eeec..f93b3e2 100644 --- a/.gitignore +++ b/.gitignore @@ -29,5 +29,5 @@ CMakeUserPresets.json sliced.vtk contour.vtk -testData/ +testData/cases/*/*.vtk build diff --git a/CMakePresets.json b/CMakePresets.json index 93003ba..bb86429 100644 --- a/CMakePresets.json +++ b/CMakePresets.json @@ -10,6 +10,7 @@ "type": "FILEPATH", "value": "$env{VCPKG_ROOT}/scripts/buildsystems/vcpkg.cmake" }, + "CMAKE_BUILD_TYPE": "Release", "TESSELLATOR_ENABLE_CGAL": false } }, @@ -26,13 +27,14 @@ { "name": "gnu", "displayName": "GNU g++ compiler", + "binaryDir": "build-rls/", "inherits": "default" }, { "name": "gnu-cgal", "displayName": "GNU g++ compiler with CGAL", "inherits": "gnu", - "binaryDir": "build-cgal/", + "binaryDir": "build-rls-cgal/", "cacheVariables": { "TESSELLATOR_ENABLE_CGAL": true } diff --git a/README.md b/README.md index 7b1e9a2..97f884c 100644 --- a/README.md +++ b/README.md @@ -95,11 +95,12 @@ This contains the information about the mesh file(s). You can specify a single o - `objects`: An array of object definitions. Each object can have: - `filename`: (required) The mesh file name, relative to the JSON file location - `group`: (optional) Group name for the object (defaults to filename without extension) + - `ghost`: (optional boolean, default: false) Excludes the object from cross-object decisions while still meshing and exporting it normally - `mesher`: (optional) Override the global mesher settings for this specific object ```json "objects": [ - {"filename": "object1.stl", "group": "group1"}, + {"filename": "object1.stl", "group": "group1", "ghost": true}, {"filename": "object2.stl", "group": "group2", "mesher": {"type": "conformal"}} ] ``` @@ -118,8 +119,14 @@ For **staircase** mesher: - `splitHexahedra`: (boolean, default: false) Splits filled volumes into one conforming hexahedron per occupied grid cell For **conformal** mesher: -- `edgePoints`: Controls edge point snapping behavior -- `forbiddenLength`: Minimum length threshold for snapping +- `edgePoints`: (non-negative integer, default: `0`) Number of evenly spaced + candidate snap points added along each grid edge. These points are placed in + the portion of the edge outside the endpoint exclusion regions defined by + `forbiddenLength`. Set to `0` to add no interior edge points. +- `forbiddenLength`: (number, default: `0.0`) Fraction of each grid edge kept + clear next to both endpoints when placing or snapping to edge points. It must + not exceed `0.5`. +- `staircaseSharedCells`: (boolean, default: true) Selectively staircases cells occupied by this conformal object and another object **Global options:** - `exportGrid`: (boolean, default: true) Controls whether to export the grid file @@ -135,12 +142,26 @@ Example with staircase mesher and compression enabled: } ``` +### `` +This optional entry controls how multi-object results are written: + +- `singleFile`: (boolean, default: false) Writes all objects to one + `{basename}.tessellator.vtk` file instead of separate per-object mesh files. + Object groups remain identifiable through the `group` and `groupNames` cell + attributes. Group names must be unique. + +```json + "output": { + "singleFile": true + } +``` + Example with conformal mesher: ```json "mesher": { "type": "conformal", "options": { - "edgePoints": true, + "edgePoints": 3, "forbiddenLength": 0.001 } } @@ -150,6 +171,7 @@ Example with conformal mesher: The tessellator generates output files with the following naming convention: - `{group_name}.tessellator.str.vtk` - Staircase meshed object - `{group_name}.tessellator.cmsh.vtk` - Conformal meshed object +- `{basename}.tessellator.vtk` - Combined multi-object mesh when `output.singleFile` is true - `{basename}.tessellator.grid.vtk` - Grid file (if `exportGrid` is true) ### Complete Example diff --git a/src/app/launcher.cpp b/src/app/launcher.cpp index 5221e49..e878249 100644 --- a/src/app/launcher.cpp +++ b/src/app/launcher.cpp @@ -4,8 +4,10 @@ #include "meshers/MesherBase.h" #include "meshers/StaircaseMesher.h" #include "meshers/ConformalMesher.h" +#include "core/Staircaser.h" #include "utils/GridTools.h" #include "utils/MeshTools.h" +#include "utils/RedundancyCleaner.h" #include #include @@ -14,8 +16,10 @@ #include #include #include +#include #include #include +#include namespace meshlib::app { @@ -59,6 +63,7 @@ std::vector readObjectsFromJSON(const nlohmann::json& fileData if (obj.contains("volume")){ objDef.isVolume = obj["volume"]; } + objDef.ghost = obj.value("ghost", false); if (obj.contains("mesher")) { objDef.mesherOverride = obj["mesher"]; } @@ -71,6 +76,7 @@ std::vector readObjectsFromJSON(const nlohmann::json& fileData if (fileData["object"].contains("volume")){ objDef.isVolume = fileData["object"]["volume"]; } + objDef.ghost = fileData["object"].value("ghost", false); if (fileData.contains("mesher")) { objDef.mesherOverride = fileData["mesher"]; } @@ -82,6 +88,14 @@ std::vector readObjectsFromJSON(const nlohmann::json& fileData return objects; } +bool readSingleFileOutputOption(const nlohmann::json& fileData) +{ + if (!fileData.contains("output")) { + return false; + } + return fileData["output"].value("singleFile", false); +} + Mesh readMesh(const nlohmann::json& fileData, const std::filesystem::path& folderPath, const ObjectDefinition& objDef) { std::filesystem::path meshObjectPath = folderPath / objDef.filename; @@ -179,8 +193,13 @@ meshlib::meshers::ConformalMesherOptions readConformalMesherOptions(const nlohma res.volumeGroups.insert(0); } if (mesherConfig.contains("options")) { - res.snapperOptions.edgePoints = mesherConfig["options"]["edgePoints"]; - res.snapperOptions.forbiddenLength = mesherConfig["options"]["forbiddenLength"]; + const auto& options = mesherConfig["options"]; + res.snapperOptions.edgePoints = options.value( + "edgePoints", res.snapperOptions.edgePoints); + res.snapperOptions.forbiddenLength = options.value( + "forbiddenLength", res.snapperOptions.forbiddenLength); + res.staircaseSharedCells = options.value( + "staircaseSharedCells", res.staircaseSharedCells); } return res; } @@ -218,6 +237,105 @@ std::unique_ptr buildMesher(const Mesh& in, const } } +namespace { + +struct MeshedObject { + ObjectDefinition definition; + std::string mesherType; + std::string extension; + bool staircaseSharedCells = false; + Mesh mesh; +}; + +Mesh mergeObjectMeshes(const std::vector& objects) +{ + Mesh combined; + if (objects.empty()) { + return combined; + } + + combined.grid = objects.front().mesh.grid; + for (const auto& object : objects) { + if (object.mesh.grid != combined.grid) { + throw std::runtime_error("Cannot combine object meshes with different grids."); + } + if (object.mesh.groups.size() != 1) { + throw std::runtime_error( + "Each object must produce exactly one mesh group for combined processing."); + } + + Mesh namedMesh = object.mesh; + namedMesh.groups.front().name = object.definition.group; + if (combined.groups.empty()) { + combined = std::move(namedMesh); + } else { + utils::meshTools::mergeMeshAsNewGroup(combined, namedMesh); + } + } + utils::RedundancyCleaner::fuseCoords(combined); + utils::RedundancyCleaner::cleanCoords(combined); + return combined; +} + +void validateCombinedGroupNames(const std::vector& objects) +{ + std::set groupNames; + for (const auto& object : objects) { + if (!groupNames.insert(object.group).second) { + throw std::runtime_error( + "Combined output requires unique object group names; duplicate group: " + + object.group); + } + } +} + +void staircaseSharedConformalCells(std::vector& objects) +{ + const bool hasEnabledConformalObject = std::any_of( + objects.begin(), objects.end(), [](const MeshedObject& object) { + return !object.definition.ghost && + object.mesherType == conformal_mesher && + object.staircaseSharedCells; + }); + const auto participatingObjectCount = std::count_if( + objects.begin(), objects.end(), [](const MeshedObject& object) { + return !object.definition.ghost; + }); + if (!hasEnabledConformalObject || participatingObjectCount < 2) { + return; + } + + Mesh relativeCombined = mergeObjectMeshes(objects); + utils::meshTools::convertToRelativeCoordinates(relativeCombined); + std::set ghostGroups; + for (GroupId groupId = 0; groupId < objects.size(); ++groupId) { + if (objects[groupId].definition.ghost) { + ghostGroups.insert(groupId); + } + } + const auto sharedCells = meshlib::meshers::ConformalMesher::cellsSharedByGroups( + relativeCombined, ghostGroups); + if (sharedCells.empty()) { + return; + } + + for (auto& object : objects) { + if (object.definition.ghost || object.mesherType != conformal_mesher || + !object.staircaseSharedCells) { + continue; + } + + Mesh relativeMesh = object.mesh; + utils::meshTools::convertToRelativeCoordinates(relativeMesh); + object.mesh = meshlib::core::Staircaser{relativeMesh}.getSelectiveMesh( + sharedCells, meshlib::core::Staircaser::GapsFillingType::Insert); + object.mesh.groups.front().name = object.definition.group; + utils::meshTools::convertToAbsoluteCoordinates(object.mesh); + } +} + +} // namespace + int launcher(int argc, const char* argv[]) { po::options_description desc("Allowed options"); @@ -247,9 +365,13 @@ int launcher(int argc, const char* argv[]) std::vector objects = readObjectsFromJSON(inputFileData); std::filesystem::path outputFolder = getFolder(inputFileName); auto basename = getBasename(inputFileName); + const bool singleFileOutput = readSingleFileOutputOption(inputFileData); + if (singleFileOutput) { + validateCombinedGroupNames(objects); + } - Mesh firstMesh; - bool first = true; + std::vector meshedObjects; + meshedObjects.reserve(objects.size()); for (const auto& objDef : objects) { std::cout << "\n-- Processing object: " << objDef.filename << " (group: " << objDef.group << ")" << std::endl; @@ -258,20 +380,46 @@ int launcher(int argc, const char* argv[]) auto mesher = buildMesher(mesh, inputFileData, objDef); Mesh resultMesh = mesher->mesh(); + if (resultMesh.groups.size() == 1) { + resultMesh.groups.front().name = objDef.group; + } - if (first) { - firstMesh = resultMesh; - first = false; + const auto mesherType = readMesherType(inputFileData, objDef.mesherOverride); + bool staircaseSharedCells = false; + if (mesherType == conformal_mesher) { + staircaseSharedCells = readConformalMesherOptions( + inputFileData, objDef.isVolume, objDef.mesherOverride) + .staircaseSharedCells; } + meshedObjects.push_back({ + objDef, + mesherType, + readExtension(inputFileData, objDef.mesherOverride), + staircaseSharedCells, + std::move(resultMesh) + }); + } + + staircaseSharedConformalCells(meshedObjects); - auto extension = readExtension(inputFileData, objDef.mesherOverride); - std::string outputFileName = objDef.group + ".tessellator." + extension + ".vtk"; - exportMeshToVTU(outputFolder / outputFileName, resultMesh); + if (singleFileOutput && !meshedObjects.empty()) { + const auto outputFileName = basename + ".tessellator.vtk"; + exportMeshToVTU( + outputFolder / outputFileName, mergeObjectMeshes(meshedObjects)); std::cout << "-- Exported: " << outputFileName << std::endl; + } else { + for (const auto& object : meshedObjects) { + const std::string outputFileName = object.definition.group + + ".tessellator." + object.extension + ".vtk"; + exportMeshToVTU(outputFolder / outputFileName, object.mesh); + std::cout << "-- Exported: " << outputFileName << std::endl; + } } - if (!first && readExportGridOption(inputFileData, std::nullopt)) { - exportGridToVTU(outputFolder / (basename + ".tessellator.grid.vtk"), firstMesh.grid); + if (!meshedObjects.empty() && readExportGridOption(inputFileData, std::nullopt)) { + exportGridToVTU( + outputFolder / (basename + ".tessellator.grid.vtk"), + meshedObjects.front().mesh.grid); std::cout << "-- Exported grid: " << basename << ".tessellator.grid.vtk" << std::endl; } diff --git a/src/app/launcher.h b/src/app/launcher.h index 48823c3..c6f98e6 100644 --- a/src/app/launcher.h +++ b/src/app/launcher.h @@ -16,13 +16,15 @@ struct ObjectDefinition { std::string filename; std::string group; bool isVolume = false; + bool ghost = false; std::optional mesherOverride; }; int launcher(int argc, const char* argv[]); Grid parseGridFromJSON(const nlohmann::json& fileData); std::vector readObjectsFromJSON(const nlohmann::json& fileData); +bool readSingleFileOutputOption(const nlohmann::json& fileData); Mesh readMesh(const std::string& fn, const ObjectDefinition& objDef); std::unique_ptr buildMesher(const Mesh& in, const nlohmann::json& fileData, const ObjectDefinition& objDef); -} \ No newline at end of file +} diff --git a/src/app/vtkIO.cpp b/src/app/vtkIO.cpp index 71f00b5..d9c887c 100644 --- a/src/app/vtkIO.cpp +++ b/src/app/vtkIO.cpp @@ -24,6 +24,19 @@ const char* GROUPS_TAG_NAME = "group"; namespace meshlib::vtkIO { +namespace +{ + +std::string normalizeLegacyVTKString(std::string value) +{ + if (!value.empty() && value.back() == '\r') { + value.pop_back(); + } + return value; +} + +} + vtkSmartPointer vtkPolyDataToVTU(vtkPolyData* polyData) { vtkNew appendFilter; @@ -148,11 +161,17 @@ Mesh vtuToMesh(vtkUnstructuredGrid* vtu) if (vtu->GetCellData()->HasArray(GROUPS_TAG_NAME)) { vtkIntArray* groupsDataArray = vtkIntArray::SafeDownCast(vtu->GetCellData()->GetArray(GROUPS_TAG_NAME)); + vtkStringArray* groupNamesDataArray = vtkStringArray::SafeDownCast( + vtu->GetCellData()->GetAbstractArray("groupNames")); mesh.groups.resize(groupsDataArray->GetRange()[1] + 1); for (vtkIdType i = 0; i < vtu->GetNumberOfCells(); i++) { auto g = groupsDataArray->GetValue(i); mesh.groups[g].elements.push_back( vtkCellToElement(vtu->GetCell(i))); + if (groupNamesDataArray != nullptr && mesh.groups[g].name.empty()) { + mesh.groups[g].name = normalizeLegacyVTKString( + groupNamesDataArray->GetValue(i)); + } } } else { mesh.groups.resize(1); diff --git a/src/core/Smoother.cpp b/src/core/Smoother.cpp index 9dc38ea..663fb38 100644 --- a/src/core/Smoother.cpp +++ b/src/core/Smoother.cpp @@ -8,6 +8,7 @@ #include #include +#include #include #include @@ -167,5 +168,45 @@ Smoother::Smoother(const Mesh& mesh, const SmootherOptions& opts) : meshTools::checkNoCellsAreCrossed(mesh_); } +Mesh Smoother::retriangulatePlanarPatches( + const Mesh& mesh, + double featureDetectionAngle) +{ + Mesh result = mesh; + SmootherTools tools(result.grid); + for (auto& group : result.groups) { + std::vector patches; + for (const auto& cell : tools.buildCellElemMap(group.elements, result.coordinates)) { + ElementsView triangles; + std::copy_if( + cell.second.begin(), cell.second.end(), + std::back_inserter(triangles), + [&](const Element* element) { + return element->isTriangle() + && !Geometry::isDegenerate( + Geometry::asTriV(*element, result.coordinates)); + }); + if (triangles.empty()) { + continue; + } + for (const auto& patch : Geometry::buildDisjointSmoothSets( + triangles, result.coordinates, featureDetectionAngle)) { + if (CoordGraph(patch).getBoundAndInteriorVertices().second.empty()) { + patches.push_back(patch); + } + } + } + for (const auto& patch : patches) { + tools.retriangulatePlanarPatch( + group.elements, result.coordinates, patch); + } + } + + RedundancyCleaner::removeDegenerateElements(result); + RedundancyCleaner::cleanCoords(result); + meshTools::checkNoCellsAreCrossed(result); + return result; +} + } } diff --git a/src/core/Smoother.h b/src/core/Smoother.h index 7bc2bb4..168e1c7 100644 --- a/src/core/Smoother.h +++ b/src/core/Smoother.h @@ -17,6 +17,10 @@ class Smoother { Smoother(const Mesh&, const SmootherOptions& opts = SmootherOptions()); Mesh getMesh() const { return mesh_; } + static Mesh retriangulatePlanarPatches( + const Mesh& mesh, + double featureDetectionAngle = SmootherOptions().featureDetectionAngle); + private: SmootherOptions opts_; SmootherTools sT_; @@ -29,4 +33,4 @@ class Smoother { }; } -} \ No newline at end of file +} diff --git a/src/core/SmootherTools.cpp b/src/core/SmootherTools.cpp index dc46199..c0c1e3b 100644 --- a/src/core/SmootherTools.cpp +++ b/src/core/SmootherTools.cpp @@ -1,12 +1,15 @@ #include "Smoother.h" +#include "Slicer.h" #include "utils/CoordGraph.h" #include "utils/ElemGraph.h" +#include "utils/ConvexHull.h" #include "utils/RedundancyCleaner.h" #include "utils/Geometry.h" #include "utils/Tools.h" #include +#include #include #include #include @@ -411,6 +414,69 @@ void SmootherTools::remeshBoundary( } } +void SmootherTools::retriangulatePlanarPatch( + Elements& es, + const Coordinates& cs, + const ElementsView& patch) +{ + if (patch.empty()) { + return; + } + + const CoordGraph graph(patch); + const auto [boundary, interior] = graph.getBoundAndInteriorVertices(); + if (!interior.empty() || !patchIsPlanar(cs, patch)) { + return; + } + + const auto cycles = graph.getBoundaryGraph().findCycles(); + if (cycles.size() != 1) { + return; + } + + const auto referenceNormal = Geometry::normal( + Geometry::asTriV(*patch[0], cs)); + const auto path = ConvexHull(&cs).get(boundary, referenceNormal); + if (path.size() != boundary.size()) { + return; + } + + using Edge = std::pair; + std::set boundaryEdges; + std::set hullEdges; + for (std::size_t i = 0; i < path.size(); ++i) { + hullEdges.insert(std::minmax(path[i], path[(i + 1) % path.size()])); + } + for (const auto& line : graph.getBoundaryGraph().getEdgesAsLines()) { + boundaryEdges.insert(std::minmax(line.vertices[0], line.vertices[1])); + } + if (hullEdges != boundaryEdges) { + return; + } + + auto remeshedElements = Slicer::buildTrianglesFromPath(cs, path); + if (remeshedElements.size() != patch.size() + || std::any_of( + remeshedElements.begin(), remeshedElements.end(), + [&](const Element& element) { + return Geometry::isDegenerate(Geometry::asTriV(element, cs)); + })) { + return; + } + + for (auto& element : remeshedElements) { + if (hasWrongOrientation(*patch[0], element, cs)) { + reorientSingleElement(element); + } + } + + const std::lock_guard lock(writingElements_); + for (std::size_t i = 0; i < patch.size(); ++i) { + const ElementId elementId = patch[i] - &es.front(); + es[elementId] = remeshedElements[i]; + } +} + void SmootherTools::remeshWithNoInteriorPoints( Elements& es, const Coordinates& cs, diff --git a/src/core/SmootherTools.h b/src/core/SmootherTools.h index 3ab957b..6235ffd 100644 --- a/src/core/SmootherTools.h +++ b/src/core/SmootherTools.h @@ -75,6 +75,11 @@ class SmootherTools : public utils::GridTools { const Coordinates& meshCs, const ElementsView& patch); + void retriangulatePlanarPatch( + Elements& es, + const Coordinates& cs, + const ElementsView& patch); + SingularIds buildSingularIds( const Elements& es, const Coordinates& cs, diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index f51f2af..93f01c4 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -17,6 +17,9 @@ Staircaser::Staircaser(const Mesh& inputMesh) : GridTools(inputMesh.grid) mesh_.coordinates.reserve(inputMesh.coordinates.size() * 2); mesh_.groups.resize(inputMesh.groups.size()); + for (std::size_t groupId = 0; groupId < inputMesh.groups.size(); ++groupId) { + mesh_.groups[groupId].name = inputMesh.groups[groupId].name; + } } @@ -155,23 +158,33 @@ Mesh Staircaser::getSelectiveMesh(const std::set& cellsToStructure, GapsFi auto cellElemMap = buildCellElemMap(inputGroup.elements, inputMesh_.coordinates); + CoordinateMap coordinateMap = buildCoordinateMap(mesh_.coordinates); for (const auto& c : cellsToStructure) { if (!cellElemMap.count(c)) { continue; } for (const auto e: cellElemMap.at(c)) { - if (e->isLine()) { + if (e->isNode()) { + this->processNodeAndAddToGroup( + *e, inputMesh_.coordinates, mesh_.coordinates, meshGroup); + } + else if (e->isLine()) { this->processLineAndAddToGroup(*e, inputMesh_.coordinates, mesh_.coordinates, meshGroup); } else if (e->isTriangle()) { this->processTriangleAndAddToGroup(*e, inputMesh_.coordinates, meshGroup); } + else { + auto [newElement, ignoredBoundaryCoordinates] = + obtainNewIndexForElement(*e, {}, coordinateMap); + meshGroup.elements.push_back(std::move(newElement)); + } } } meshGroup.elements.reserve(meshGroup.elements.size() + inputGroup.elements.size()); std::map> cellElemMap_withNewElements; - CoordinateMap coordinateMap = buildCoordinateMap(mesh_.coordinates); + coordinateMap = buildCoordinateMap(mesh_.coordinates); for (const auto& [cell, elements] : cellElemMap) { if (cellsToStructure.count(cell) ) { diff --git a/src/meshers/ConformalMesher.cpp b/src/meshers/ConformalMesher.cpp index dca455c..9d8b790 100644 --- a/src/meshers/ConformalMesher.cpp +++ b/src/meshers/ConformalMesher.cpp @@ -137,6 +137,34 @@ std::set ConformalMesher::cellsWithAVertexInAnEdgeForbiddenRegion(const Me return res; } +std::set ConformalMesher::cellsSharedByGroups( + const Mesh& mesh, const std::set& ignoredGroups) +{ + const GridTools gridTools(mesh.grid); + std::map> groupsByCell; + + for (GroupId groupId = 0; groupId < mesh.groups.size(); ++groupId) { + if (ignoredGroups.count(groupId) != 0) { + continue; + } + const auto elementsByCell = gridTools.buildCellElemMap( + mesh.groups[groupId].elements, mesh.coordinates); + for (const auto& [cell, elements] : elementsByCell) { + if (!elements.empty()) { + groupsByCell[cell].insert(groupId); + } + } + } + + std::set sharedCells; + for (const auto& [cell, groups] : groupsByCell) { + if (groups.size() > 1) { + sharedCells.insert(cell); + } + } + return sharedCells; +} + std::set mergeCellSets(const std::set& a, const std::set& b) { std::set res; @@ -188,6 +216,10 @@ Mesh ConformalMesher::mesh() const log("Snapping.", 1); res = Snapper(res, opts_.snapperOptions).getMesh(); + + log("Retriangulating planar patches.", 1); + res = Smoother::retriangulatePlanarPatches( + res, smootherOpts.featureDetectionAngle); logNumberOfTriangles(countMeshElementsIf(res, isTriangle)); // Find cells which break conformal FDTD rules. diff --git a/src/meshers/ConformalMesher.h b/src/meshers/ConformalMesher.h index 26a9858..96c6b52 100644 --- a/src/meshers/ConformalMesher.h +++ b/src/meshers/ConformalMesher.h @@ -23,6 +23,8 @@ class ConformalMesher : public MesherBase { static std::set findNonConformalCells(const Mesh& mesh); static std::set cellsWithMoreThanAVertexInsideEdge(const Mesh& mesh); static std::set cellsWithMoreThanAPathPerFace(const Mesh& mesh); + static std::set cellsSharedByGroups( + const Mesh& mesh, const std::set& ignoredGroups = {}); static std::set cellsWithInteriorDisconnectedPatches(const Mesh& mesh); static std::set cellsWithAVertexInAnEdgeForbiddenRegion(const Mesh& mesh); private: diff --git a/src/meshers/ConformalMesherOptions.h b/src/meshers/ConformalMesherOptions.h index 5a38c18..1680146 100644 --- a/src/meshers/ConformalMesherOptions.h +++ b/src/meshers/ConformalMesherOptions.h @@ -9,7 +9,7 @@ namespace meshlib::meshers { class ConformalMesherOptions : public MesherBaseOptions { public: core::SnapperOptions snapperOptions; - // std::set volumeGroups{}; + bool staircaseSharedCells = true; }; } diff --git a/src/meshers/OffgridMesher.cpp b/src/meshers/OffgridMesher.cpp index e7c8914..278756b 100644 --- a/src/meshers/OffgridMesher.cpp +++ b/src/meshers/OffgridMesher.cpp @@ -58,6 +58,7 @@ void OffgridMesher::process(Mesh& mesh) const if (opts_.snap) { log("Snapping.", 1); mesh = Snapper(mesh, opts_.snapperOptions).getMesh(); + mesh = Smoother::retriangulatePlanarPatches(mesh); logNumberOfTriangles(countMeshElementsIf(mesh, isTriangle)); } } diff --git a/src/utils/MeshTools.cpp b/src/utils/MeshTools.cpp index 00e0f1b..f1b7856 100644 --- a/src/utils/MeshTools.cpp +++ b/src/utils/MeshTools.cpp @@ -419,7 +419,7 @@ void mergeMeshAsNewGroup(Mesh& lMesh, const Mesh& iMesh) lMesh.coordinates.insert(lMesh.coordinates.end(), iMesh.coordinates.begin(), iMesh.coordinates.end()); - lMesh.groups.push_back(Group()); + lMesh.groups.push_back(Group(iMesh.groups.front().name, {})); mergeGroup(lMesh.groups.back(), iMesh.groups.front(), coordCount); } @@ -490,4 +490,4 @@ Mesh extractGroupsByName(const Mesh& mesh, const std::vector& group return result; } -} \ No newline at end of file +} diff --git a/test/app/launcherTest.cpp b/test/app/launcherTest.cpp index 2cf01ab..7ee2139 100644 --- a/test/app/launcherTest.cpp +++ b/test/app/launcherTest.cpp @@ -3,6 +3,7 @@ #include "app/launcher.h" #include "types/Mesh.h" #include "meshers/StaircaseMesher.h" +#include "meshers/ConformalMesher.h" #include #include @@ -99,6 +100,57 @@ TEST_F(LauncherTest, buildsStaircasedMesherWithSplitHexahedra) EXPECT_TRUE(staircase.getOptions().splitHexahedra); } +TEST_F(LauncherTest, singleFileOutputIsDisabledByDefault) +{ + EXPECT_FALSE(readSingleFileOutputOption(nlohmann::json::object())); + EXPECT_TRUE(readSingleFileOutputOption({ + {"output", {{"singleFile", true}}} + })); +} + +TEST_F(LauncherTest, conformalMesherStaircasesSharedCellsByDefault) +{ + meshlib::Mesh meshMock; + meshMock.grid = { + std::vector{0, 1}, + std::vector{0, 1}, + std::vector{0, 1} + }; + const nlohmann::json config = { + {"mesher", {{"type", "conformal"}}} + }; + + ObjectDefinition object; + auto mesher = buildMesher(meshMock, config, object); + const auto& conformal = + dynamic_cast(*mesher); + + EXPECT_TRUE(conformal.getOptions().staircaseSharedCells); +} + +TEST_F(LauncherTest, conformalSharedCellStaircasingCanBeDisabled) +{ + meshlib::Mesh meshMock; + meshMock.grid = { + std::vector{0, 1}, + std::vector{0, 1}, + std::vector{0, 1} + }; + const nlohmann::json config = { + {"mesher", { + {"type", "conformal"}, + {"options", {{"staircaseSharedCells", false}}} + }} + }; + + ObjectDefinition object; + auto mesher = buildMesher(meshMock, config, object); + const auto& conformal = + dynamic_cast(*mesher); + + EXPECT_FALSE(conformal.getOptions().staircaseSharedCells); +} + TEST_F(LauncherTest, builds_staircased_mesher_without_compression) { meshlib::Mesh meshMock; @@ -251,10 +303,12 @@ TEST_F(LauncherTest, readObjectsFromJSON_basic) EXPECT_EQ(objects[0].filename, "sphere.stl"); EXPECT_EQ(objects[0].group, "sphere_group"); EXPECT_TRUE(objects[0].isVolume); + EXPECT_FALSE(objects[0].ghost); EXPECT_FALSE(objects[0].mesherOverride.has_value()); EXPECT_EQ(objects[1].filename, "cone.stl"); EXPECT_EQ(objects[1].group, "cone_group"); EXPECT_FALSE(objects[1].isVolume); + EXPECT_FALSE(objects[1].ghost); EXPECT_FALSE(objects[1].mesherOverride.has_value()); } @@ -317,6 +371,51 @@ TEST_F(LauncherTest, readObjectsFromJSON_legacyFormat) EXPECT_TRUE(objects[0].isVolume); } +TEST_F(LauncherTest, readObjectsFromJSON_solenoid) +{ + std::ifstream input("testData/cases/solenoid/solenoid.tessellator.json"); + nlohmann::json config; + input >> config; + + const auto objects = readObjectsFromJSON(config); + + ASSERT_EQ(objects.size(), 5); + EXPECT_EQ(objects[0].filename, "solenoid.vtu"); + EXPECT_EQ(objects[0].group, "Solenoid"); + EXPECT_FALSE(objects[0].ghost); + ASSERT_TRUE(objects[0].mesherOverride.has_value()); + EXPECT_EQ(objects[0].mesherOverride.value()["type"], "conformal"); + EXPECT_EQ(objects[1].filename, "BC.vtu"); + EXPECT_EQ(objects[1].group, "BC"); + EXPECT_TRUE(objects[1].ghost); + EXPECT_EQ(objects[2].filename, "Generator.vtu"); + EXPECT_EQ(objects[2].group, "Generator"); + EXPECT_FALSE(objects[2].ghost); + EXPECT_EQ(objects[3].filename, "wire_left.vtu"); + EXPECT_EQ(objects[3].group, "Wire_left"); + EXPECT_FALSE(objects[3].ghost); + EXPECT_EQ(objects[4].filename, "wire_right.vtu"); + EXPECT_EQ(objects[4].group, "Wire_right"); + EXPECT_FALSE(objects[4].ghost); + + meshlib::Mesh meshMock; + meshMock.grid = { + std::vector{0, 1}, + std::vector{0, 1}, + std::vector{0, 1} + }; + auto solenoidMesher = buildMesher(meshMock, config, objects[0]); + EXPECT_NE( + dynamic_cast(solenoidMesher.get()), + nullptr); + for (std::size_t objectIndex = 1; objectIndex < objects.size(); ++objectIndex) { + auto mesher = buildMesher(meshMock, config, objects[objectIndex]); + EXPECT_NE( + dynamic_cast(mesher.get()), + nullptr); + } +} + TEST_F(LauncherTest, builds_staircased_mesher_with_override) { meshlib::Mesh meshMock; @@ -352,6 +451,92 @@ TEST_F(LauncherTest, launches_multiObject_basic) EXPECT_EQ(exitCode, EXIT_SUCCESS); } +TEST_F(LauncherTest, launchesMultiObjectIntoSingleGroupedFile) +{ + const auto temp = std::filesystem::temp_directory_path(); + const auto input = temp / "tessellator_single_file.tessellator.json"; + const auto output = temp / "tessellator_single_file.tessellator.vtk"; + const auto gridOutput = temp / "tessellator_single_file.tessellator.grid.vtk"; + const auto sphereOutput = temp / "sphere_group.tessellator.str.vtk"; + const auto coneOutput = temp / "cone_group.tessellator.str.vtk"; + std::filesystem::remove(output); + std::filesystem::remove(gridOutput); + std::filesystem::remove(sphereOutput); + std::filesystem::remove(coneOutput); + const nlohmann::json config = { + {"grid", { + {"numberOfCells", {10, 10, 10}}, + {"boundingBox", {{-100, -100, -100}, {100, 100, 100}}} + }}, + {"mesher", {{"type", "staircase"}}}, + {"output", {{"singleFile", true}}}, + {"objects", { + { + {"filename", std::filesystem::absolute( + "testData/cases/multiObject/sphere.stl").string()}, + {"volume", true}, + {"group", "sphere_group"} + }, + { + {"filename", std::filesystem::absolute( + "testData/cases/multiObject/cone.stl").string()}, + {"group", "cone_group"} + } + }} + }; + { + std::ofstream stream(input); + stream << config; + } + + const std::string inputString = input.string(); + const char* av[] = {nullptr, "-i", inputString.c_str()}; + EXPECT_EQ(launcher(3, av), EXIT_SUCCESS); + + ASSERT_TRUE(std::filesystem::exists(output)); + EXPECT_FALSE(std::filesystem::exists(sphereOutput)); + EXPECT_FALSE(std::filesystem::exists(coneOutput)); + { + std::ifstream stream(output); + const std::string contents{ + std::istreambuf_iterator(stream), std::istreambuf_iterator()}; + EXPECT_NE(contents.find("groupNames"), std::string::npos); + EXPECT_NE(contents.find("sphere_group"), std::string::npos); + EXPECT_NE(contents.find("cone_group"), std::string::npos); + } + + std::filesystem::remove(output); + std::filesystem::remove(gridOutput); + std::filesystem::remove(input); +} + +TEST_F(LauncherTest, rejectsDuplicateGroupNamesInSingleFileOutput) +{ + const auto input = std::filesystem::temp_directory_path() + / "tessellator_duplicate_groups.json"; + const nlohmann::json config = { + {"grid", { + {"numberOfCells", {1, 1, 1}}, + {"boundingBox", {{0, 0, 0}, {1, 1, 1}}} + }}, + {"output", {{"singleFile", true}}}, + {"objects", { + {{"filename", "first.stl"}, {"group", "duplicate"}}, + {{"filename", "second.stl"}, {"group", "duplicate"}} + }} + }; + { + std::ofstream stream(input); + stream << config; + } + const std::string inputString = input.string(); + const char* av[] = {nullptr, "-i", inputString.c_str()}; + + EXPECT_THROW(launcher(3, av), std::runtime_error); + + std::filesystem::remove(input); +} + TEST_F(LauncherTest, launches_multiObject_mixedMesher) { int ac = 3; @@ -378,3 +563,12 @@ TEST_F(LauncherTest, launches_multiObject_sameFileMultipleGroups) EXPECT_NO_THROW(exitCode = launcher(ac, av)); EXPECT_EQ(exitCode, EXIT_SUCCESS); } + +TEST_F(LauncherTest, launches_solenoid_multiObject_case) +{ + int ac = 3; + const char* av[] = { NULL, "-i", "testData/cases/solenoid/solenoid.tessellator.json"}; + int exitCode; + EXPECT_NO_THROW(exitCode = launcher(ac, av)); + EXPECT_EQ(exitCode, EXIT_SUCCESS); +} diff --git a/test/app/vtkIOTest.cpp b/test/app/vtkIOTest.cpp index 4e89dc5..482f5d4 100644 --- a/test/app/vtkIOTest.cpp +++ b/test/app/vtkIOTest.cpp @@ -3,6 +3,9 @@ #include "app/vtkIO.h" #include "utils/GridTools.h" +#include +#include + using namespace meshlib::vtkIO; class VTKIOTest : public ::testing::Test @@ -43,6 +46,49 @@ TEST_F(VTKIOTest, readElementTypes) EXPECT_TRUE(m.groups[0].elements[2].isTriangle()); } +TEST_F(VTKIOTest, readsSolenoidBoundaryConditionAsRectangle) +{ + const auto mesh = readInputMesh("testData/cases/solenoid/BC.vtu"); + + ASSERT_EQ(mesh.coordinates.size(), 4); + ASSERT_EQ(mesh.groups.size(), 1); + ASSERT_EQ(mesh.groups[0].elements.size(), 2); + EXPECT_TRUE(mesh.groups[0].elements[0].isTriangle()); + EXPECT_TRUE(mesh.groups[0].elements[1].isTriangle()); + EXPECT_DOUBLE_EQ(mesh.coordinates[0][0], -5.8284271247); + EXPECT_DOUBLE_EQ(mesh.coordinates[0][1], -5.7573593129); + EXPECT_DOUBLE_EQ(mesh.coordinates[0][2], 10.0); + EXPECT_DOUBLE_EQ(mesh.coordinates[2][0], 19.627416998); + EXPECT_DOUBLE_EQ(mesh.coordinates[2][1], 19.69848481); + EXPECT_DOUBLE_EQ(mesh.coordinates[2][2], 10.0); +} + +TEST_F(VTKIOTest, solenoidTrianglesHaveConsistentWindingAcrossSharedEdges) +{ + const auto mesh = readInputMesh("testData/cases/solenoid/solenoid.vtu"); + ASSERT_EQ(mesh.groups.size(), 1); + + using Edge = std::pair; + std::map> edgeUses; + for (const auto& triangle : mesh.groups[0].elements) { + ASSERT_TRUE(triangle.isTriangle()); + for (std::size_t vertex = 0; vertex < triangle.vertices.size(); ++vertex) { + const Edge oriented{ + triangle.vertices[vertex], + triangle.vertices[(vertex + 1) % triangle.vertices.size()]}; + edgeUses[std::minmax(oriented.first, oriented.second)].push_back(oriented); + } + } + + for (const auto& [edge, uses] : edgeUses) { + ASSERT_LE(uses.size(), 2) << "Non-manifold edge " << edge.first << '-' << edge.second; + if (uses.size() == 2) { + EXPECT_EQ(uses[0], Edge(uses[1].second, uses[1].first)) + << "Inconsistent winding at edge " << edge.first << '-' << edge.second; + } + } +} + TEST_F(VTKIOTest, exportAndReadHexahedron) { meshlib::Mesh mesh; @@ -67,6 +113,68 @@ TEST_F(VTKIOTest, exportAndReadHexahedron) EXPECT_EQ(mesh.coordinates, result.coordinates); } +TEST_F(VTKIOTest, exportAndReadGroupNames) +{ + meshlib::Mesh mesh; + mesh.coordinates = { + meshlib::Coordinate({0.0, 0.0, 0.0}), + meshlib::Coordinate({1.0, 0.0, 0.0}) + }; + mesh.groups = { + meshlib::Group("first", { + meshlib::Element({0}, meshlib::Element::Type::Node)}), + meshlib::Group("second", { + meshlib::Element({1}, meshlib::Element::Type::Node)}) + }; + + const auto filename = std::filesystem::temp_directory_path() + / "tessellator_group_names_roundtrip.vtu"; + exportMeshToVTU(filename, mesh); + const auto result = readInputMesh(filename); + std::filesystem::remove(filename); + + ASSERT_EQ(result.groups.size(), 2); + EXPECT_EQ(result.groups[0].name, "first"); + EXPECT_EQ(result.groups[1].name, "second"); +} + +TEST_F(VTKIOTest, readsGroupNamesFromLegacyVTKWithCRLFLineEndings) +{ + const auto filename = std::filesystem::temp_directory_path() + / "tessellator_crlf_group_names.vtu"; + { + std::ofstream stream(filename, std::ios::binary); + stream + << "# vtk DataFile Version 3.0\r\n" + << "group names with CRLF\r\n" + << "ASCII\r\n" + << "DATASET UNSTRUCTURED_GRID\r\n" + << "POINTS 2 float\r\n" + << "0 0 0\r\n" + << "1 0 0\r\n" + << "CELLS 2 4\r\n" + << "1 0\r\n" + << "1 1\r\n" + << "CELL_TYPES 2\r\n" + << "1\r\n" + << "1\r\n" + << "CELL_DATA 2\r\n" + << "FIELD FieldData 2\r\n" + << "group 1 2 int\r\n" + << "0 1\r\n" + << "groupNames 1 2 string\r\n" + << "first\r\n" + << "second\r\n"; + } + + const auto result = readInputMesh(filename); + std::filesystem::remove(filename); + + ASSERT_EQ(result.groups.size(), 2); + EXPECT_EQ(result.groups[0].name, "first"); + EXPECT_EQ(result.groups[1].name, "second"); +} + TEST_F(VTKIOTest, exportGridToVTU) { meshlib::Grid grid; diff --git a/test/core/SmootherTest.cpp b/test/core/SmootherTest.cpp index fce3dac..6d39bb9 100644 --- a/test/core/SmootherTest.cpp +++ b/test/core/SmootherTest.cpp @@ -7,6 +7,9 @@ #include "utils/MeshTools.h" #include "core/Slicer.h" +#include +#include + #if APP_LOADED #include "app/vtkIO.h" #endif @@ -90,6 +93,203 @@ TEST_F(SmootherTest, touching_by_single_point) EXPECT_EQ(1, countMeshElementsIf(r, isTriangle)); } +TEST_F(SmootherTest, retriangulatesMirroredPlanarPatchesWithMirroredDiagonals) +{ + Mesh mesh; + mesh.grid = utils::GridTools::buildCartesianGrid(0.0, 1.0, 2); + mesh.coordinates = { + Relative({0.0, 0.0, 0.0}), + Relative({0.0, 0.0, 1.0}), + Relative({1.0, 1.0, 1.0}), + Relative({1.0, 1.0, 0.0}) + }; + mesh.groups = {Group(), Group()}; + mesh.groups[0].elements = { + Element({0, 3, 1}, Element::Type::Surface), + Element({1, 3, 2}, Element::Type::Surface) + }; + mesh.groups[1].elements = { + Element({0, 2, 3}, Element::Type::Surface), + Element({0, 1, 2}, Element::Type::Surface) + }; + + const auto result = Smoother::retriangulatePlanarPatches(mesh); + const auto internalEdge = [](const Elements& elements) { + std::map, std::size_t> edgeUses; + for (const auto& triangle : elements) { + for (std::size_t vertex = 0; vertex < triangle.vertices.size(); ++vertex) { + ++edgeUses[std::minmax( + triangle.vertices[vertex], + triangle.vertices[(vertex + 1) % triangle.vertices.size()])]; + } + } + return std::find_if( + edgeUses.begin(), edgeUses.end(), + [](const auto& edgeUse) { return edgeUse.second == 2; })->first; + }; + + EXPECT_EQ((std::pair{0, 2}), + internalEdge(result.groups[0].elements)); + EXPECT_EQ((std::pair{1, 3}), + internalEdge(result.groups[1].elements)); + for (const auto& triangle : result.groups[0].elements) { + EXPECT_GT(Geometry::normal(Geometry::asTriV(triangle, result.coordinates))[X], 0.0); + } + for (const auto& triangle : result.groups[1].elements) { + EXPECT_LT(Geometry::normal(Geometry::asTriV(triangle, result.coordinates))[X], 0.0); + } +} + +TEST_F(SmootherTest, retriangulatesTrianglesAndPreservesNonTriangularElements) +{ + Mesh mesh; + mesh.grid = utils::GridTools::buildCartesianGrid(0.0, 1.0, 2); + mesh.coordinates = { + Relative({0.0, 0.0, 0.0}), + Relative({0.0, 0.0, 1.0}), + Relative({1.0, 1.0, 1.0}), + Relative({1.0, 1.0, 0.0}) + }; + mesh.groups = {Group()}; + mesh.groups[0].elements = { + Element({0}, Element::Type::Node), + Element({0, 1}, Element::Type::Line), + Element({0, 1, 2, 3}, Element::Type::Surface), + Element({0, 3, 1}, Element::Type::Surface), + Element({1, 3, 2}, Element::Type::Surface) + }; + + const auto result = Smoother::retriangulatePlanarPatches(mesh); + + ASSERT_EQ(result.groups[0].elements.size(), 5); + EXPECT_EQ(result.groups[0].elements[0], mesh.groups[0].elements[0]); + EXPECT_EQ(result.groups[0].elements[1], mesh.groups[0].elements[1]); + EXPECT_EQ(result.groups[0].elements[2], mesh.groups[0].elements[2]); + std::map, std::size_t> triangleEdgeUses; + for (const auto& element : result.groups[0].elements) { + if (!element.isTriangle()) { + continue; + } + for (std::size_t vertex = 0; vertex < element.vertices.size(); ++vertex) { + ++triangleEdgeUses[std::minmax( + element.vertices[vertex], + element.vertices[(vertex + 1) % element.vertices.size()])]; + } + } + EXPECT_EQ((triangleEdgeUses[std::pair{0, 2}]), 2); + EXPECT_EQ((triangleEdgeUses[std::pair{1, 3}]), 0); +} + +TEST_F(SmootherTest, preservesConcavePlanarPatch) +{ + Mesh mesh; + mesh.grid = utils::GridTools::buildCartesianGrid(0.0, 1.0, 2); + mesh.coordinates = { + Relative({0.1, 0.1, 0.5}), Relative({0.9, 0.1, 0.5}), + Relative({0.9, 0.5, 0.5}), Relative({0.5, 0.5, 0.5}), + Relative({0.5, 0.9, 0.5}), Relative({0.1, 0.9, 0.5}) + }; + mesh.groups = {Group()}; + mesh.groups[0].elements = { + Element({0, 1, 3}), Element({1, 2, 3}), + Element({0, 3, 5}), Element({3, 4, 5}) + }; + + const auto result = Smoother::retriangulatePlanarPatches(mesh); + + EXPECT_EQ(result.groups[0].elements, mesh.groups[0].elements); +} + +TEST_F(SmootherTest, preservesDegenerateTriangles) +{ + Mesh mesh; + mesh.grid = utils::GridTools::buildCartesianGrid(0.0, 1.0, 2); + mesh.coordinates = { + Relative({0.1, 0.1, 0.5}), + Relative({0.5, 0.5, 0.5}), + Relative({0.9, 0.9, 0.5}) + }; + mesh.groups = {Group()}; + mesh.groups[0].elements = {Element({0, 1, 2})}; + + const auto result = Smoother::retriangulatePlanarPatches(mesh); + + ASSERT_EQ(result.groups[0].elements.size(), 1); + EXPECT_EQ(result.groups[0].elements[0], mesh.groups[0].elements[0]); +} + +TEST_F(SmootherTest, preservesPlanarPatchWithMultipleBoundaryCycles) +{ + Mesh mesh; + mesh.grid = utils::GridTools::buildCartesianGrid(0.0, 1.0, 2); + mesh.coordinates = { + Relative({0.1, 0.1, 0.5}), Relative({0.9, 0.1, 0.5}), + Relative({0.9, 0.9, 0.5}), Relative({0.1, 0.9, 0.5}), + Relative({0.35, 0.35, 0.5}), Relative({0.65, 0.35, 0.5}), + Relative({0.65, 0.65, 0.5}), Relative({0.35, 0.65, 0.5}) + }; + mesh.groups = {Group()}; + mesh.groups[0].elements = { + Element({0, 1, 5}), Element({0, 5, 4}), + Element({1, 2, 6}), Element({1, 6, 5}), + Element({2, 3, 7}), Element({2, 7, 6}), + Element({3, 0, 4}), Element({3, 4, 7}) + }; + + const auto result = Smoother::retriangulatePlanarPatches(mesh); + + EXPECT_EQ(result.groups[0].elements, mesh.groups[0].elements); +} + +TEST_F(SmootherTest, retriangulatesPlanarPatchesInAdjacentCellsIndependently) +{ + Mesh mesh; + mesh.grid = utils::GridTools::buildCartesianGrid(0.0, 2.0, 3); + mesh.coordinates = { + Relative({0.1, 0.1, 0.5}), Relative({0.1, 0.9, 0.5}), + Relative({1.0, 0.9, 0.5}), Relative({1.0, 0.1, 0.5}), + Relative({1.9, 0.9, 0.5}), Relative({1.9, 0.1, 0.5}) + }; + mesh.groups = {Group()}; + mesh.groups[0].elements = { + Element({0, 3, 1}), Element({1, 3, 2}), + Element({3, 5, 2}), Element({2, 5, 4}) + }; + + const auto result = Smoother::retriangulatePlanarPatches(mesh); + std::set> edges; + for (const auto& triangle : result.groups[0].elements) { + for (std::size_t vertex = 0; vertex < triangle.vertices.size(); ++vertex) { + edges.insert(std::minmax( + triangle.vertices[vertex], + triangle.vertices[(vertex + 1) % triangle.vertices.size()])); + } + } + + EXPECT_EQ(edges.count({0, 2}), 1); + EXPECT_EQ(edges.count({3, 4}), 1); + EXPECT_EQ(edges.count({1, 3}), 0); + EXPECT_EQ(edges.count({2, 5}), 0); + EXPECT_NO_THROW(meshTools::checkNoCellsAreCrossed(result)); +} + +TEST_F(SmootherTest, rejectsElementsCrossingCellBoundaries) +{ + Mesh mesh; + mesh.grid = utils::GridTools::buildCartesianGrid(0.0, 2.0, 3); + mesh.coordinates = { + Relative({0.5, 0.2, 0.5}), + Relative({1.5, 0.2, 0.5}), + Relative({0.5, 0.8, 0.5}) + }; + mesh.groups = {Group()}; + mesh.groups[0].elements = {Element({0, 1, 2})}; + + EXPECT_THROW( + Smoother::retriangulatePlanarPatches(mesh), + std::runtime_error); +} + #if APP_LOADED TEST_F(SmootherTest, preserves_topological_closedness_for_alhambra) @@ -152,4 +352,4 @@ TEST_F(SmootherTest, preserves_topological_closedness_for_sphere) #endif -} \ No newline at end of file +} diff --git a/test/core/StaircaserTest.cpp b/test/core/StaircaserTest.cpp index baa63b1..57c8b5f 100644 --- a/test/core/StaircaserTest.cpp +++ b/test/core/StaircaserTest.cpp @@ -2842,3 +2842,25 @@ TEST_F(StaircaserTest, selectiveStructurer_SplitLinesWithNeighborTriangle) } } + +TEST_F(StaircaserTest, selectiveStaircasingPreservesStructuredElementsAndGroupName) +{ + Mesh mesh; + mesh.grid = GridTools::buildCartesianGrid(0.0, 1.0, 2); + mesh.coordinates = { + Relative({0.1, 0.1, 0.5}), + Relative({0.9, 0.1, 0.5}), + Relative({0.9, 0.9, 0.5}), + Relative({0.1, 0.9, 0.5}) + }; + mesh.groups = { + Group("structured", {Element({0, 1, 2, 3}, Element::Type::Surface)}) + }; + + const auto result = Staircaser{mesh}.getSelectiveMesh({Cell({0, 0, 0})}); + + ASSERT_EQ(result.groups.size(), 1); + EXPECT_EQ(result.groups.front().name, "structured"); + ASSERT_EQ(result.groups.front().elements.size(), 1); + EXPECT_TRUE(result.groups.front().elements.front().isQuad()); +} diff --git a/test/meshers/ConformalMesherTest.cpp b/test/meshers/ConformalMesherTest.cpp index b13db5d..736ee8a 100644 --- a/test/meshers/ConformalMesherTest.cpp +++ b/test/meshers/ConformalMesherTest.cpp @@ -42,6 +42,66 @@ class ConformalMesherTest : public ::testing::Test { #endif }; +TEST_F(ConformalMesherTest, findsCellOccupiedByDifferentGroups) +{ + Mesh mesh; + mesh.grid = buildUnitLengthGrid(0.5); + mesh.coordinates = { + Relative({0.25, 0.25, 0.25}), + Relative({0.75, 0.75, 0.75}), + Relative({1.25, 0.25, 0.25}) + }; + mesh.groups = { + Group("first", {Element({0}, Element::Type::Node)}), + Group("second", {Element({1}, Element::Type::Node)}), + Group("third", {Element({2}, Element::Type::Node)}) + }; + + const auto result = ConformalMesher::cellsSharedByGroups(mesh); + + EXPECT_EQ(result, std::set({Cell({0, 0, 0})})); +} + +TEST_F(ConformalMesherTest, includesBothCellsWhenGroupsMeetOnCellFace) +{ + Mesh mesh; + mesh.grid = buildUnitLengthGrid(0.5); + mesh.coordinates = { + Relative({1.0, 0.25, 0.25}), + Relative({1.0, 0.75, 0.75}) + }; + mesh.groups = { + Group("first", {Element({0}, Element::Type::Node)}), + Group("second", {Element({1}, Element::Type::Node)}) + }; + + const auto result = ConformalMesher::cellsSharedByGroups(mesh); + + EXPECT_EQ(result, std::set({Cell({0, 0, 0}), Cell({1, 0, 0})})); +} + +TEST_F(ConformalMesherTest, ignoresSelectedGroupsWhenFindingSharedCells) +{ + Mesh mesh; + mesh.grid = buildUnitLengthGrid(0.5); + mesh.coordinates = { + Relative({0.25, 0.25, 0.25}), + Relative({0.75, 0.75, 0.75}), + Relative({1.25, 0.25, 0.25}), + Relative({1.75, 0.75, 0.75}) + }; + mesh.groups = { + Group("first", {Element({0}, Element::Type::Node)}), + Group("ghost", {Element({1}, Element::Type::Node)}), + Group("second", {Element({2}, Element::Type::Node)}), + Group("third", {Element({3}, Element::Type::Node)}) + }; + + const auto result = ConformalMesher::cellsSharedByGroups(mesh, {1}); + + EXPECT_EQ(result, std::set({Cell({1, 0, 0})})); +} + TEST_F(ConformalMesherTest, cellsWithMoreThanAVertexPerEdge_1) { // This is non-conformal. @@ -490,4 +550,4 @@ TEST_F(ConformalMesherTest, thinCylinder) -} \ No newline at end of file +} diff --git a/testData/cases/solenoid/BC.vtu b/testData/cases/solenoid/BC.vtu new file mode 100644 index 0000000..b2cbcda --- /dev/null +++ b/testData/cases/solenoid/BC.vtu @@ -0,0 +1,17 @@ +# vtk DataFile Version 5.1 +vtk output +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 4 double +-5.8284271247 -5.7573593129 10 +19.627416998 -5.7573593129 10 +19.627416998 19.69848481 10 +-5.8284271247 19.69848481 10 +CELLS 3 6 +OFFSETS vtktypeint64 +0 3 6 +CONNECTIVITY vtktypeint64 +0 1 2 0 2 3 +CELL_TYPES 2 +5 +5 diff --git a/testData/cases/solenoid/Generator.vtu b/testData/cases/solenoid/Generator.vtu new file mode 100644 index 0000000..c8d6a03 --- /dev/null +++ b/testData/cases/solenoid/Generator.vtu @@ -0,0 +1,13 @@ +# vtk DataFile Version 5.1 +vtk output +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 1 double +-1.3322676296e-15 14.142135624 0 +CELLS 2 1 +OFFSETS vtktypeint64 +0 1 +CONNECTIVITY vtktypeint64 +0 +CELL_TYPES 1 +1 diff --git a/testData/cases/solenoid/solenoid.tessellator.json b/testData/cases/solenoid/solenoid.tessellator.json new file mode 100644 index 0000000..1acd561 --- /dev/null +++ b/testData/cases/solenoid/solenoid.tessellator.json @@ -0,0 +1,26 @@ +{ + "_version": "0.1", + "_format": "Tessellator Data File in JSON format", + "output": { + "singleFile": true + }, + "grid": { + "numberOfCells": [74, 75, 60], + "boundingBox": [ + [-34.141999999999996, -25.756999999999998, -20.0], + [39.858, 49.243, 40.0] + ] + }, + "objects": [ + {"filename": "solenoid.vtu", "group": "Solenoid", "mesher": {"type": "conformal"}}, + {"filename": "BC.vtu", "group": "BC", "ghost": true}, + {"filename": "Generator.vtu", "group": "Generator"}, + {"filename": "wire_left.vtu", "group": "Wire_left"}, + {"filename": "wire_right.vtu", "group": "Wire_right"} + ], + "mesher": { + "options": { + "exportGrid": true + } + } +} diff --git a/testData/cases/solenoid/solenoid.vtu b/testData/cases/solenoid/solenoid.vtu new file mode 100644 index 0000000..f514d30 --- /dev/null +++ b/testData/cases/solenoid/solenoid.vtu @@ -0,0 +1,40 @@ +# vtk DataFile Version 5.1 +vtk output +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 16 double +-1.9539925233e-14 28.28427124 20 -14.14213562 14.14213562 20 1.4142135144e-15 -1.4142135144e-15 20 +14.14213562 14.14213562 20 0.70710676908 13.43502903 0 -6.3639612198 6.3639612198 0 +0 0 0 7.0710678101 7.0710678101 0 14.14213562 14.14213562 0 +-14.14213562 14.14213562 0 -7.0710678101 21.21320343 0 -1.9539925233e-14 28.28427124 0 +-0.70710676908 14.84924221 0 6.3639612198 21.920310974 0 7.7781744003 20.506095886 0 +-7.7781744003 7.7781744003 0 +CELLS 17 48 +OFFSETS vtktypeint64 +0 3 6 9 12 15 18 21 24 +27 30 33 36 39 42 45 48 +CONNECTIVITY vtktypeint64 +0 1 2 0 2 3 4 6 5 +4 7 6 2 6 7 3 7 8 +3 2 7 1 10 9 0 11 10 +0 10 1 11 12 10 11 13 12 +14 7 4 14 8 7 10 15 9 +10 12 15 +CELL_TYPES 16 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 +5 + diff --git a/testData/cases/solenoid/wire_left.vtu b/testData/cases/solenoid/wire_left.vtu new file mode 100644 index 0000000..033599f --- /dev/null +++ b/testData/cases/solenoid/wire_left.vtu @@ -0,0 +1,13 @@ +# vtk DataFile Version 5.1 +vtk output +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 2 double +-0.7071 14.8492 0 -0 14.1421 0 +CELLS 2 2 +OFFSETS vtktypeint64 +0 2 +CONNECTIVITY vtktypeint64 +0 1 +CELL_TYPES 1 +3 diff --git a/testData/cases/solenoid/wire_right.vtu b/testData/cases/solenoid/wire_right.vtu new file mode 100644 index 0000000..c9f062e --- /dev/null +++ b/testData/cases/solenoid/wire_right.vtu @@ -0,0 +1,13 @@ +# vtk DataFile Version 5.1 +vtk output +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 2 double +-0 14.1421 0 0.7071 13.435 0 +CELLS 2 2 +OFFSETS vtktypeint64 +0 2 +CONNECTIVITY vtktypeint64 +0 1 +CELL_TYPES 1 +3