diff --git a/.github/workflows/build-and-test.yml b/.github/workflows/build-and-test.yml index 6f78829..c6b8acd 100644 --- a/.github/workflows/build-and-test.yml +++ b/.github/workflows/build-and-test.yml @@ -19,9 +19,10 @@ jobs: matrix: preset: [ {"os": windows-2022, "name": "msbuild"}, - {"os": ubuntu-latest, "name": "gnu"} + {"os": ubuntu-latest, "name": "gnu"}, + {"os": ubuntu-latest, "name": "gnu-cgal"} ] - build-type: ["Debug", "Release"] + build-type: ["Release"] fail-fast: false @@ -85,7 +86,7 @@ jobs: run: build/bin/${{matrix.build-type}}/tessellator_tests.exe - name: Ubuntu configure and build - if: matrix.preset.name=='gnu' + if: matrix.preset.os=='ubuntu-latest' run: | configured=0 for attempt in 1 2 3; do @@ -99,5 +100,5 @@ jobs: cmake --build build -j - name: Ubuntu Run tests - if: matrix.preset.name=='gnu' + if: matrix.preset.os=='ubuntu-latest' run: build/bin/tessellator_tests diff --git a/.gitignore b/.gitignore index 1fd2b13..130eeec 100644 --- a/.gitignore +++ b/.gitignore @@ -14,6 +14,7 @@ src/*.json .vs/ .vscode/settings.json +.vscode/settings.dev.json .vscode/launch.json .vscode/tasks.json diff --git a/.vscode/settings.dev.json b/.vscode/settings.dev.json index 5303634..9d190ee 100644 --- a/.vscode/settings.dev.json +++ b/.vscode/settings.dev.json @@ -8,12 +8,8 @@ "C_Cpp.intelliSenseEngine": "disabled", "clangd.arguments": ["-log=verbose", "-pretty", - "--background-index", + "--background-index", "--query-driver=/usr/bin/g++" ], - "clangd.fallbackFlags": ["-std=c++17"], - "lldb.showDisassembly": "never", - "lldb.launch.preRunCommands": [ - "settings set target.import-std-module true" - ] + "clangd.fallbackFlags": ["-std=c++17"] } diff --git a/CLAUDE.md b/CLAUDE.md index fa6ab49..f9e2bac 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -64,6 +64,10 @@ Uses **CMake 3.20+** with presets and vcpkg for dependency management. cmake --preset gnu -S . -B build cmake --build build -j +# Optional CGAL algorithms +cmake --preset gnu-cgal +cmake --build --preset gnu-cgal -j + # Windows (MSBuild) cmake --preset msbuild -S . -B build cmake --build build --config Release -j @@ -72,7 +76,7 @@ cmake --build build --config Release -j ### CMake Options - `TESSELLATOR_ENABLE_TESTS` (ON by default) – Build test suite -- `TESSELLATOR_ENABLE_CGAL` (ON by default) – Enable CGAL-based geometry operations +- `TESSELLATOR_ENABLE_CGAL` (OFF by default) – Enable CGAL-based geometry operations - `TESSELLATOR_EXECUTION_POLICIES` (OFF by default) – Parallel execution policies ### Dependencies @@ -94,7 +98,7 @@ To set up locally, create a `CMakeUserPreset.json` file: "VCPKG_ROOT": "~/workspace/vcpkg/" }, "cacheVariables": { - "TESSELLATOR_ENABLE_CGAL": true + "TESSELLATOR_ENABLE_CGAL": false }, "inherits": "gnu" } diff --git a/CMakeLists.txt b/CMakeLists.txt index 577007a..5d3a1b0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -27,7 +27,7 @@ set(CMAKE_EXPORT_COMPILE_COMMANDS OFF) option(TESSELLATOR_ENABLE_TESTS "Compile tests" ON) -option(TESSELLATOR_ENABLE_CGAL "Compile using CGAL library" ON) +option(TESSELLATOR_ENABLE_CGAL "Compile using CGAL library" OFF) option(TESSELLATOR_EXECUTION_POLICIES OFF) option(DOCKER_EXPORT_COMPILE_COMMANDS "Information for the 'clangd' linter extension" OFF) option(TESSELLATOR_LOAD_APP "Compile app" OFF) diff --git a/CMakePresets.json b/CMakePresets.json index 769d22e..93003ba 100644 --- a/CMakePresets.json +++ b/CMakePresets.json @@ -10,7 +10,7 @@ "type": "FILEPATH", "value": "$env{VCPKG_ROOT}/scripts/buildsystems/vcpkg.cmake" }, - "TESSELLATOR_ENABLE_CGAL": true + "TESSELLATOR_ENABLE_CGAL": false } }, { @@ -26,9 +26,17 @@ { "name": "gnu", "displayName": "GNU g++ compiler", - "generator": "Ninja", "inherits": "default" }, + { + "name": "gnu-cgal", + "displayName": "GNU g++ compiler with CGAL", + "inherits": "gnu", + "binaryDir": "build-cgal/", + "cacheVariables": { + "TESSELLATOR_ENABLE_CGAL": true + } + }, { "name": "gnu-dbg", "displayName": "GNU g++ compiler - Debug", @@ -49,7 +57,7 @@ "CMAKE_PREFIX_PATH": "/usr/local", "CMAKE_FIND_ROOT_PATH": "/usr/local", "TESSELLATOR_ENABLE_TESTS": "ON", - "TESSELLATOR_ENABLE_CGAL": "ON", + "TESSELLATOR_ENABLE_CGAL": "OFF", "DOCKER_EXPORT_COMPILE_COMMANDS": "ON" }, "environment": { @@ -68,7 +76,7 @@ "CMAKE_PREFIX_PATH": "/usr/local", "CMAKE_FIND_ROOT_PATH": "/usr/local", "TESSELLATOR_ENABLE_TESTS": "ON", - "TESSELLATOR_ENABLE_CGAL": "ON", + "TESSELLATOR_ENABLE_CGAL": "OFF", "DOCKER_EXPORT_COMPILE_COMMANDS": "ON", "CMAKE_CXX_FLAGS_INIT": "-g3" }, @@ -89,6 +97,10 @@ { "name": "docker-dbg", "configurePreset": "docker-dbg" + }, + { + "name": "gnu-cgal", + "configurePreset": "gnu-cgal" } ] -} \ No newline at end of file +} diff --git a/README.md b/README.md index c920511..7b1e9a2 100644 --- a/README.md +++ b/README.md @@ -18,6 +18,7 @@ Tessellator is a mesher focused on generate meshes and data structures which are ## Compilation When using presets, make sure to define the environment variable `VCPKG_ROOT` to your `vcpkg` installation. +The standard presets build without CGAL; use the `gnu-cgal` preset for the optional CGAL algorithms. This can be done using a `CMakeUserPreset.json` file, for example: ```json @@ -33,7 +34,7 @@ This can be done using a `CMakeUserPreset.json` file, for example: "VCPKG_ROOT": "~/workspace/vcpkg/" }, "cacheVariables": { - "TESSELLATOR_ENABLE_CGAL": true + "TESSELLATOR_ENABLE_CGAL": false }, "inherits": "gnu" } @@ -114,6 +115,7 @@ This optional entry configures the meshing algorithm and its options. If not spe For **staircase** mesher: - `compress`: (boolean, default: false) Enables surface compression to merge adjacent coplanar quads into larger surfaces +- `splitHexahedra`: (boolean, default: false) Splits filled volumes into one conforming hexahedron per occupied grid cell For **conformal** mesher: - `edgePoints`: Controls edge point snapping behavior diff --git a/src/app/launcher.cpp b/src/app/launcher.cpp index 6ce1de5..5221e49 100644 --- a/src/app/launcher.cpp +++ b/src/app/launcher.cpp @@ -150,17 +150,22 @@ meshlib::meshers::StaircaseMesherOptions readStaircaseMesherOptions(const nlohma } meshlib::meshers::StaircaseMesherOptions res; - - res.isVolume = isVolume; + if (isVolume){ + res.volumeGroups.insert(0); + } if (mesherConfig.contains("options") && mesherConfig["options"].contains("compress")) { res.compress = mesherConfig["options"]["compress"]; } + if (mesherConfig.contains("options") && + mesherConfig["options"].contains("splitHexahedra")) { + res.splitHexahedra = mesherConfig["options"]["splitHexahedra"]; + } return res; } -meshlib::meshers::ConformalMesherOptions readConformalMesherOptions(const nlohmann::json& fileData, const std::optional& override) +meshlib::meshers::ConformalMesherOptions readConformalMesherOptions(const nlohmann::json& fileData, bool isVolume, const std::optional& override) { nlohmann::json mesherConfig; if (override.has_value()) { @@ -170,6 +175,9 @@ meshlib::meshers::ConformalMesherOptions readConformalMesherOptions(const nlohma } meshlib::meshers::ConformalMesherOptions res; + if (isVolume){ + res.volumeGroups.insert(0); + } if (mesherConfig.contains("options")) { res.snapperOptions.edgePoints = mesherConfig["options"]["edgePoints"]; res.snapperOptions.forbiddenLength = mesherConfig["options"]["forbiddenLength"]; @@ -204,7 +212,7 @@ std::unique_ptr buildMesher(const Mesh& in, const readStaircaseMesherOptions(fileData, objDef.isVolume, objDef.mesherOverride) }); } else if (mesherType == meshlib::app::conformal_mesher) { - return std::make_unique(meshlib::meshers::ConformalMesher{in, readConformalMesherOptions(fileData, objDef.mesherOverride)}); + return std::make_unique(meshlib::meshers::ConformalMesher{in, readConformalMesherOptions(fileData, objDef.isVolume, objDef.mesherOverride)}); } else { throw std::runtime_error("Unsupported mesher type"); } diff --git a/src/app/vtkIO.cpp b/src/app/vtkIO.cpp index 0bdc367..71f00b5 100644 --- a/src/app/vtkIO.cpp +++ b/src/app/vtkIO.cpp @@ -3,7 +3,10 @@ #include #include #include +#include +#include #include +#include #include #include #include @@ -76,6 +79,8 @@ Element vtkCellToElement(vtkCell* cell) vtkVertex* vertex = nullptr; vtkLine* line = nullptr; vtkTriangle* triangle = nullptr; + vtkTetra* tetra = nullptr; + vtkHexahedron* hexahedron = nullptr; switch (cell->GetCellType()) { case VTK_VERTEX: @@ -102,6 +107,26 @@ Element vtkCellToElement(vtkCell* cell) }; elem.type = meshlib::Element::Type::Surface; break; + + case VTK_TETRA: + tetra = vtkTetra::SafeDownCast(cell); + elem.vertices = { + CoordinateId(tetra->GetPointIds()->GetId(0)), + CoordinateId(tetra->GetPointIds()->GetId(1)), + CoordinateId(tetra->GetPointIds()->GetId(2)), + CoordinateId(tetra->GetPointIds()->GetId(3)) + }; + elem.type = meshlib::Element::Type::Volume; + break; + + case VTK_HEXAHEDRON: + hexahedron = vtkHexahedron::SafeDownCast(cell); + elem.vertices.reserve(8); + for (vtkIdType id = 0; id < 8; ++id) { + elem.vertices.push_back(CoordinateId(hexahedron->GetPointIds()->GetId(id))); + } + elem.type = meshlib::Element::Type::Volume; + break; } return elem; @@ -131,6 +156,7 @@ Mesh vtuToMesh(vtkUnstructuredGrid* vtu) } } else { mesh.groups.resize(1); + auto k = vtu->GetNumberOfCells(); mesh.groups[0].elements.reserve(vtu->GetNumberOfCells()); for (vtkIdType i = 0; i < vtu->GetNumberOfCells(); i++) { mesh.groups[0].elements.push_back( @@ -205,6 +231,12 @@ vtkSmartPointer elementsToVTU(const Mesh& mesh) } else if (elem.isNode()) { cellTypes.push_back(VTK_VERTEX); cell = vtkSmartPointer::New(); + } else if (elem.isTetrahedron()) { + cellTypes.push_back(VTK_TETRA); + cell = vtkSmartPointer::New(); + } else if (elem.isHexahedron()) { + cellTypes.push_back(VTK_HEXAHEDRON); + cell = vtkSmartPointer::New(); } else { throw std::runtime_error("Unsupported element type"); } diff --git a/src/cgal/filler/Filler.cpp b/src/cgal/filler/Filler.cpp index ae72d0a..60c23bf 100644 --- a/src/cgal/filler/Filler.cpp +++ b/src/cgal/filler/Filler.cpp @@ -253,6 +253,29 @@ void sliceAlignedByGrid( ); } +void sliceAlignedByGridAndRemove( + Filler::GridSlices& slices, + const Polyhedron& m, + const Grid& g, + const Priority& priority) +{ + + auto polygons{ buildGridPlanesPolygons(makeFacesCCWOriented(m), g)}; + const std::array axis{ X, Y, Z }; + + std::for_each( +#ifdef TESSELLATOR_EXECUTION_POLICIES + std::execution::par, +#endif + axis.begin(), axis.end(), + [&](const auto& x) { + for (const auto& [i, polygon] : polygons[x]) { + slices[x][i].remove(polygon, priority); + } + } + ); +} + Priority Filler::getGroupPriority(const GroupId& gId) const { if (gId < groupPriorities_.size()) { @@ -508,7 +531,8 @@ FillerPolyhedrons buildFillerPolyhedrons( Filler::Filler( const Mesh& volumeMesh, const Mesh& surfaceMesh, - const std::vector& groupPriorities) + const std::vector& groupPriorities, + const FillerMode& fillerMode) { utils::meshTools::checkNoNullAreasExist(volumeMesh); utils::meshTools::checkNoNullAreasExist(surfaceMesh); @@ -547,7 +571,11 @@ Filler::Filler( log("Slicing surfaces", 2); sliceNonAlignedByGrid(slices_, fP.surfaces, grid_, pr, SlicingMode::Surface); log("Slicing aligned", 2); - sliceAlignedByGrid(slices_, fP.aligned, grid_, pr); + if (fillerMode == FillerMode::insideAndOutside) { + sliceAlignedByGrid(slices_, fP.aligned, grid_, pr); + } else if (fillerMode == FillerMode::onlyInside){ + sliceAlignedByGridAndRemove(slices_, fP.aligned, grid_, pr); + } log("Building segments arrays", 2); buildSegmentsArray(segmentsArray_, fP.aligned, grid_, pr); buildSegmentsArray(segmentsArray_, fP.volumes, grid_, pr); diff --git a/src/cgal/filler/Filler.h b/src/cgal/filler/Filler.h index 2b2d5af..ed87558 100644 --- a/src/cgal/filler/Filler.h +++ b/src/cgal/filler/Filler.h @@ -9,17 +9,26 @@ namespace meshlib::cgal::filler { +enum class FillerMode{ + insideAndOutside, + onlyInside +}; + + class Filler { -public: + public: using Slices = std::map; using GridSlices = std::array; using SegmentsArray = std::map; using GridSegmentsArray = std::array; - + + FillerMode mode = FillerMode::insideAndOutside; + Filler( const Mesh& volumeMesh, const Mesh& surfaceMesh = Mesh(), - const std::vector& groupPriorities = std::vector()); + const std::vector& groupPriorities = std::vector(), + const FillerMode& mode = FillerMode::insideAndOutside); Filler(const Filler&) = delete; Filler(Filler&&) = default; Filler& operator=(const Filler&) = delete; @@ -32,7 +41,7 @@ class Filler { FillingState getFillingState(const CellIndex&) const; Mesh getMeshFilling() const; - + GridSlices getSlices() const {return slices_;}; private: GridSlices slices_; GridSegmentsArray segmentsArray_; @@ -42,6 +51,7 @@ class Filler { void mergeGroupsWithSamePriority(Groups& vGroups, Groups& sGroups); + }; diff --git a/src/cgal/filler/Slice.cpp b/src/cgal/filler/Slice.cpp index 9a674d7..5133ba6 100644 --- a/src/cgal/filler/Slice.cpp +++ b/src/cgal/filler/Slice.cpp @@ -183,6 +183,21 @@ void Slice::add(const Polylines2& polylines, const Priority& pr) } } +void Slice::remove(const Polylines2& polylines, const Priority& pr) +{ + SliceData& sd = data_[pr]; + + for (const auto& p : polylines) { + if (p.size() == 1) { + continue; + } + auto r{ removeSegmentsContainedInAnyAxis(p) }; + for (const auto& rr : r){ + sd.lines.erase(std::find(sd.lines.begin(), sd.lines.end(), rr)); + } + } +} + FillingState::FillingState(const FillingType& t) : type{ t }, priority_{ 0 } @@ -261,6 +276,18 @@ void Slice::add(const HPolygonSet& polygons, const Priority& pr) removeInSuperiorPriorities(pr); } +void Slice::remove(const HPolygonSet& polygons, const Priority& pr) +{ + if (polygons.isEmpty()) { + return; + } + + SliceData& sd = data_[pr]; + sd.surfaces.difference(polygons); + + removeInSuperiorPriorities(pr); +} + void Slice::mergeLines(const Slice& lhs) { for (const auto& [pr, sd] : lhs.data_) { diff --git a/src/cgal/filler/Slice.h b/src/cgal/filler/Slice.h index b61883c..b04bf76 100644 --- a/src/cgal/filler/Slice.h +++ b/src/cgal/filler/Slice.h @@ -85,8 +85,10 @@ class Slice { FillingState getFillingState(const ArrayIndex&) const; void add(const Polylines2&, const Priority&); + void remove(const Polylines2&, const Priority&); void addAsPolygon(const Polylines2&, const Priority&); void add(const HPolygonSet&, const Priority&); + void remove(const HPolygonSet&, const Priority&); void mergeLines(const Slice& lhs); void buildSearchMap(); void buildTriangulations(); diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index faab34a..b3dc8a3 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -8,6 +8,9 @@ add_library(tessellator-core "Smoother.cpp" "SmootherTools.cpp" "Staircaser.cpp" + "VolumeFiller.cpp" + "VolumeShellExtractor.cpp" ) -target_link_libraries(tessellator-core tessellator-utils) \ No newline at end of file +target_link_libraries(tessellator-core + tessellator-utils) diff --git a/src/core/Slicer.cpp b/src/core/Slicer.cpp index fb5d7a2..8246fff 100644 --- a/src/core/Slicer.cpp +++ b/src/core/Slicer.cpp @@ -38,7 +38,6 @@ void orient(const Coordinates& coords, } } - Slicer::Slicer(const Mesh& input, const std::vector& dimensionPolicy, const SlicerOptions& opts) : GridTools(input.grid), opts_(opts) diff --git a/src/core/Slicer.h b/src/core/Slicer.h index e0bae73..b51ae15 100644 --- a/src/core/Slicer.h +++ b/src/core/Slicer.h @@ -24,8 +24,6 @@ class Slicer : public utils::GridTools { Slicer(const Mesh&, const std::vector& dimensionPolicy = {}, const SlicerOptions& opts = SlicerOptions()); Mesh getMesh() const { return mesh_; }; - - static Elements buildTrianglesFromPath(const std::vector&, const std::vector&); private: @@ -53,6 +51,11 @@ class Slicer : public utils::GridTools { const Cell&, const Coordinate&, const Cell&, const Coordinate&) const; + // IdSet buildGroupIntersectionsWithGridPlanes( + // Coordinates& sCoords, + // const std::vector& elements); + + }; } diff --git a/src/core/VolumeFiller.cpp b/src/core/VolumeFiller.cpp new file mode 100644 index 0000000..d04bbe6 --- /dev/null +++ b/src/core/VolumeFiller.cpp @@ -0,0 +1,211 @@ +#include "VolumeFiller.h" + +#include "utils/RedundancyCleaner.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace meshlib::core { + +using namespace utils; + +namespace { + +using Ray = std::array; +using Rays = std::array>, 3>; + +Cell coordinateCell(const Coordinate& coordinate, const GridTools& tools) +{ + const Relative relative = tools.getRelative(coordinate); + Cell cell; + for (Axis axis : {X, Y, Z}) { + const auto rounded = std::round(relative[axis]); + if (!GridTools::approxDir(relative[axis], rounded, 1e-7)) { + throw std::runtime_error( + "Volume shell contains a coordinate that is not on a grid vertex."); + } + cell[axis] = static_cast(rounded); + } + return cell; +} + +std::pair quadSurfel( + const Element& quad, + const Coordinates& coordinates, + const GridTools& tools) +{ + if (!quad.isQuad()) { + throw std::runtime_error( + "Volume filling requires a closed shell made exclusively of quads."); + } + + std::array cells; + for (std::size_t index = 0; index < cells.size(); ++index) { + cells[index] = coordinateCell(coordinates.at(quad.vertices[index]), tools); + } + + std::vector fixedAxes; + for (Axis axis : {X, Y, Z}) { + if (std::all_of(cells.begin(), cells.end(), [&](const Cell& cell) { + return cell[axis] == cells.front()[axis]; + })) { + fixedAxes.push_back(axis); + } + } + if (fixedAxes.size() != 1) { + throw std::runtime_error( + "Volume shell contains a quad that is not on one grid face."); + } + + const Axis axis = fixedAxes.front(); + Cell lower = cells.front(); + Cell upper = cells.front(); + for (const auto& cell : cells) { + for (Axis direction : {X, Y, Z}) { + lower[direction] = std::min(lower[direction], cell[direction]); + upper[direction] = std::max(upper[direction], cell[direction]); + } + } + for (Axis direction : {X, Y, Z}) { + if (direction != axis && upper[direction] - lower[direction] != 1) { + throw std::runtime_error( + "Volume shell contains a quad spanning more than one grid face."); + } + } + return {axis, lower}; +} + +CoordinateId findOrAddCoordinate( + Mesh& mesh, + std::map& coordinateIds, + const Cell& cell, + const GridTools& tools) +{ + const auto found = coordinateIds.find(cell); + if (found != coordinateIds.end()) { + return found->second; + } + const CoordinateId id = mesh.coordinates.size(); + mesh.coordinates.push_back(tools.getPos(GridTools::toRelative(cell))); + coordinateIds.emplace(cell, id); + return id; +} + +Element buildHexahedron( + Mesh& mesh, + std::map& coordinateIds, + const Cell& lower, + const Cell& upper, + const GridTools& tools) +{ + std::array vertices; + vertices[0] = Cell({lower[X], lower[Y], lower[Z]}); + vertices[1] = Cell({upper[X], lower[Y], lower[Z]}); + vertices[2] = Cell({upper[X], upper[Y], lower[Z]}); + vertices[3] = Cell({lower[X], upper[Y], lower[Z]}); + vertices[4] = Cell({lower[X], lower[Y], upper[Z]}); + vertices[5] = Cell({upper[X], lower[Y], upper[Z]}); + vertices[6] = Cell({upper[X], upper[Y], upper[Z]}); + vertices[7] = Cell({lower[X], upper[Y], upper[Z]}); + + Element hexahedron; + hexahedron.type = Element::Type::Volume; + for (const auto& vertex : vertices) { + hexahedron.vertices.push_back( + findOrAddCoordinate(mesh, coordinateIds, vertex, tools)); + } + return hexahedron; +} + +} + +VolumeFiller::VolumeFiller( + const Mesh& staircasedSurface, + bool splitHexahedra) : + GridTools(staircasedSurface.grid) +{ + mesh_.grid = staircasedSurface.grid; + mesh_.coordinates = staircasedSurface.coordinates; + mesh_.groups.resize(staircasedSurface.groups.size()); + std::map coordinateIds; + for (CoordinateId id = 0; id < staircasedSurface.coordinates.size(); ++id) { + coordinateIds.emplace( + coordinateCell(staircasedSurface.coordinates[id], *this), id); + } + + for (GroupId groupId = 0; groupId < staircasedSurface.groups.size(); ++groupId) { + const auto& inputGroup = staircasedSurface.groups[groupId]; + auto& outputGroup = mesh_.groups[groupId]; + outputGroup.name = inputGroup.name; + if (inputGroup.elements.empty()) { + continue; + } + + Rays rays; + for (const auto& element : inputGroup.elements) { + if (element.type != Element::Type::Surface) { + continue; + } + Axis axis; + Cell surfel; + std::tie(axis, surfel) = quadSurfel(element, staircasedSurface.coordinates, *this); + const Axis axis1 = (axis + 1) % 3; + const Axis axis2 = (axis + 2) % 3; + rays[axis][{surfel[axis1], surfel[axis2]}].insert(surfel[axis]); + } + + Axis fillAxis = X; + for (Axis axis : {Y, Z}) { + if (rays[axis].size() < rays[fillAxis].size()) { + fillAxis = axis; + } + } + const Axis axis1 = (fillAxis + 1) % 3; + const Axis axis2 = (fillAxis + 2) % 3; + for (const auto& ray : rays[fillAxis]) { + const auto& crossings = ray.second; + if (crossings.size() % 2 != 0) { + std::stringstream message; + message << "Volume shell has an odd number of crossings on a grid ray in group " + << groupId << "."; + throw std::runtime_error(message.str()); + } + for (auto crossing = crossings.begin(); crossing != crossings.end();) { + const CellDir begin = *crossing++; + const CellDir end = *crossing++; + if (begin == end) { + continue; + } + const CellDir lastBegin = splitHexahedra ? end - 1 : begin; + for (CellDir cellBegin = begin; cellBegin <= lastBegin; ++cellBegin) { + Cell lower; + Cell upper; + lower[fillAxis] = cellBegin; + upper[fillAxis] = splitHexahedra ? cellBegin + 1 : end; + lower[axis1] = ray.first[0]; + upper[axis1] = ray.first[0] + 1; + lower[axis2] = ray.first[1]; + upper[axis2] = ray.first[1] + 1; + outputGroup.elements.push_back( + buildHexahedron(mesh_, coordinateIds, lower, upper, *this)); + } + } + } + } + + RedundancyCleaner::cleanCoords(mesh_); +} + +Mesh VolumeFiller::getMesh() const +{ + return mesh_; +} + +} diff --git a/src/core/VolumeFiller.h b/src/core/VolumeFiller.h new file mode 100644 index 0000000..ab44a0f --- /dev/null +++ b/src/core/VolumeFiller.h @@ -0,0 +1,20 @@ +#pragma once + +#include "types/Mesh.h" +#include "utils/GridTools.h" + +namespace meshlib::core { + +class VolumeFiller : private utils::GridTools { +public: + explicit VolumeFiller( + const Mesh& staircasedSurface, + bool splitHexahedra = false); + + Mesh getMesh() const; + +private: + Mesh mesh_; +}; + +} diff --git a/src/core/VolumeShellExtractor.cpp b/src/core/VolumeShellExtractor.cpp new file mode 100644 index 0000000..5cfa698 --- /dev/null +++ b/src/core/VolumeShellExtractor.cpp @@ -0,0 +1,388 @@ +#include "VolumeShellExtractor.h" + +#include "utils/Geometry.h" +#include "utils/RedundancyCleaner.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace meshlib::core { + +namespace { + +using Face = std::array; +using Edge = std::array; + +Face faceKey(const Face& face) +{ + Face key = face; + std::sort(key.begin(), key.end()); + return key; +} + +Face asFace(const Element& element) +{ + return {element.vertices[0], element.vertices[1], element.vertices[2]}; +} + +Edge edgeKey(CoordinateId first, CoordinateId second) +{ + return first < second ? Edge{first, second} : Edge{second, first}; +} + +[[noreturn]] void fail(GroupId groupId, const std::string& reason) +{ + std::stringstream message; + message << "Invalid volume group " << groupId << ": " << reason; + throw std::runtime_error(message.str()); +} + +void validateVertexIds( + const Element& element, + const Coordinates& coordinates, + GroupId groupId) +{ + std::set unique; + for (const CoordinateId id : element.vertices) { + if (id >= coordinates.size()) { + fail(groupId, "an element references a coordinate outside the mesh."); + } + for (std::size_t axis = 0; axis < 3; ++axis) { + if (!std::isfinite(coordinates[id][axis])) { + fail(groupId, "an element references a non-finite coordinate."); + } + } + unique.insert(id); + } + if (unique.size() != element.vertices.size()) { + fail(groupId, "an element contains repeated vertices."); + } +} + +double signedTetrahedronVolume6( + const Coordinate& first, + const Coordinate& second, + const Coordinate& third, + const Coordinate& fourth) +{ + return ((second - first) ^ (third - first)) * (fourth - first); +} + +Face outwardFace( + Face face, + CoordinateId opposite, + const Coordinates& coordinates) +{ + const Coordinate& first = coordinates[face[0]]; + const Coordinate normal = + (coordinates[face[1]] - first) ^ (coordinates[face[2]] - first); + if (normal * (coordinates[opposite] - first) > 0.0) { + std::swap(face[1], face[2]); + } + return face; +} + +Elements extractTetrahedronBoundary( + const Group& group, + const Coordinates& coordinates, + GroupId groupId) +{ + std::set> tetrahedrons; + std::map> faces; + + for (const Element& element : group.elements) { + validateVertexIds(element, coordinates, groupId); + std::array tetrahedron{ + element.vertices[0], element.vertices[1], + element.vertices[2], element.vertices[3]}; + auto tetrahedronKey = tetrahedron; + std::sort(tetrahedronKey.begin(), tetrahedronKey.end()); + if (!tetrahedrons.insert(tetrahedronKey).second) { + fail(groupId, "it contains a duplicate tetrahedron."); + } + + if (std::abs(signedTetrahedronVolume6( + coordinates[tetrahedron[0]], coordinates[tetrahedron[1]], + coordinates[tetrahedron[2]], coordinates[tetrahedron[3]])) + <= utils::Geometry::NORM_TOLERANCE) { + fail(groupId, "it contains a degenerate tetrahedron."); + } + + for (std::size_t opposite = 0; opposite < tetrahedron.size(); ++opposite) { + Face face; + std::size_t faceIndex = 0; + for (std::size_t vertex = 0; vertex < tetrahedron.size(); ++vertex) { + if (vertex != opposite) { + face[faceIndex++] = tetrahedron[vertex]; + } + } + face = outwardFace(face, tetrahedron[opposite], coordinates); + auto& occurrences = faces[faceKey(face)]; + occurrences.push_back(face); + if (occurrences.size() > 2) { + fail(groupId, "a tetrahedron face is shared more than twice."); + } + } + } + + Elements boundary; + for (const auto& entry : faces) { + if (entry.second.size() == 1) { + const Face& face = entry.second.front(); + boundary.emplace_back( + CoordinateIds{face[0], face[1], face[2]}, Element::Type::Surface); + } + } + if (boundary.empty()) { + fail(groupId, "its tetrahedrons have no external boundary."); + } + return boundary; +} + +Elements extractSurfaceBoundary( + const Group& group, + const Coordinates& coordinates, + GroupId groupId) +{ + Elements boundary; + std::set faces; + boundary.reserve(group.elements.size()); + for (const Element& element : group.elements) { + validateVertexIds(element, coordinates, groupId); + const Face face = asFace(element); + if (!faces.insert(faceKey(face)).second) { + fail(groupId, "it contains a duplicate triangle."); + } + if (utils::Geometry::isDegenerate(utils::Geometry::asTriV(element, coordinates))) { + fail(groupId, "it contains a degenerate triangle."); + } + boundary.push_back(element); + } + return boundary; +} + +int edgeDirection(const Element& face, const Edge& edge) +{ + for (std::size_t vertex = 0; vertex < face.vertices.size(); ++vertex) { + const CoordinateId first = face.vertices[vertex]; + const CoordinateId second = face.vertices[(vertex + 1) % face.vertices.size()]; + if (first == edge[0] && second == edge[1]) { + return 1; + } + if (first == edge[1] && second == edge[0]) { + return -1; + } + } + throw std::logic_error("A face does not contain its indexed edge."); +} + +using EdgeFaces = std::map>; + +EdgeFaces buildEdgeFaces(const Elements& faces, GroupId groupId) +{ + EdgeFaces edgeFaces; + for (ElementId faceId = 0; faceId < faces.size(); ++faceId) { + const auto& vertices = faces[faceId].vertices; + for (std::size_t vertex = 0; vertex < vertices.size(); ++vertex) { + edgeFaces[edgeKey(vertices[vertex], vertices[(vertex + 1) % vertices.size()])] + .push_back(faceId); + } + } + for (const auto& entry : edgeFaces) { + if (entry.second.size() == 1) { + fail(groupId, "its surface is open."); + } + if (entry.second.size() != 2) { + fail(groupId, "its surface contains a non-manifold edge."); + } + } + return edgeFaces; +} + +void validateVertexFans( + const Elements& faces, + const EdgeFaces& edgeFaces, + GroupId groupId) +{ + std::map> vertexFaces; + std::map>> adjacency; + for (ElementId faceId = 0; faceId < faces.size(); ++faceId) { + for (CoordinateId vertex : faces[faceId].vertices) { + vertexFaces[vertex].insert(faceId); + } + } + for (const auto& entry : edgeFaces) { + const ElementId first = entry.second[0]; + const ElementId second = entry.second[1]; + for (CoordinateId vertex : entry.first) { + adjacency[vertex][first].insert(second); + adjacency[vertex][second].insert(first); + } + } + + for (const auto& entry : vertexFaces) { + const CoordinateId vertex = entry.first; + const auto& incident = entry.second; + std::set visited; + std::queue pending; + pending.push(*incident.begin()); + visited.insert(*incident.begin()); + while (!pending.empty()) { + const ElementId current = pending.front(); + pending.pop(); + for (ElementId next : adjacency[vertex][current]) { + if (visited.insert(next).second) { + pending.push(next); + } + } + } + if (visited.size() != incident.size()) { + fail(groupId, "its surface contains a non-manifold vertex."); + } + } +} + +std::vector> orientComponents( + Elements& faces, + const EdgeFaces& edgeFaces, + GroupId groupId) +{ + struct Neighbor { + ElementId face; + bool requiresFlip; + }; + std::vector> adjacency(faces.size()); + for (const auto& entry : edgeFaces) { + const ElementId first = entry.second[0]; + const ElementId second = entry.second[1]; + const bool requiresFlip = + edgeDirection(faces[first], entry.first) + == edgeDirection(faces[second], entry.first); + adjacency[first].push_back({second, requiresFlip}); + adjacency[second].push_back({first, requiresFlip}); + } + + std::vector flipped(faces.size(), -1); + std::vector> components; + for (ElementId seed = 0; seed < faces.size(); ++seed) { + if (flipped[seed] != -1) { + continue; + } + components.emplace_back(); + std::queue pending; + pending.push(seed); + flipped[seed] = 0; + while (!pending.empty()) { + const ElementId current = pending.front(); + pending.pop(); + components.back().push_back(current); + for (const Neighbor& neighbor : adjacency[current]) { + const int required = flipped[current] ^ neighbor.requiresFlip; + if (flipped[neighbor.face] == -1) { + flipped[neighbor.face] = required; + pending.push(neighbor.face); + } else if (flipped[neighbor.face] != required) { + fail(groupId, "its surface is not orientable."); + } + } + } + } + for (ElementId faceId = 0; faceId < faces.size(); ++faceId) { + if (flipped[faceId] == 1) { + std::swap(faces[faceId].vertices[1], faces[faceId].vertices[2]); + } + } + return components; +} + +void orientOutward( + Elements& faces, + const Coordinates& coordinates, + const std::vector>& components, + GroupId groupId) +{ + for (const auto& component : components) { + const Coordinate origin = coordinates[faces[component.front()].vertices[0]]; + double volume6 = 0.0; + for (ElementId faceId : component) { + const Face face = asFace(faces[faceId]); + volume6 += (coordinates[face[0]] - origin) + * ((coordinates[face[1]] - origin) ^ (coordinates[face[2]] - origin)); + } + if (std::abs(volume6) <= utils::Geometry::NORM_TOLERANCE) { + fail(groupId, "a closed component encloses zero volume."); + } + if (volume6 < 0.0) { + for (ElementId faceId : component) { + std::swap(faces[faceId].vertices[1], faces[faceId].vertices[2]); + } + } + } +} + +Elements extractGroupShell( + const Group& group, + const Coordinates& coordinates, + GroupId groupId) +{ + if (group.elements.empty()) { + return {}; + } + + const bool hasTriangles = std::any_of( + group.elements.begin(), group.elements.end(), + [](const Element& element) { return element.isTriangle(); }); + const bool hasTetrahedrons = std::any_of( + group.elements.begin(), group.elements.end(), + [](const Element& element) { return element.isTetrahedron(); }); + if (hasTriangles && hasTetrahedrons) { + fail(groupId, "triangles and tetrahedrons cannot be mixed."); + } + if (!hasTriangles && !hasTetrahedrons) { + fail(groupId, "only triangles or tetrahedrons are supported."); + } + if (!std::all_of( + group.elements.begin(), group.elements.end(), + [hasTriangles](const Element& element) { + return hasTriangles ? element.isTriangle() : element.isTetrahedron(); + })) { + fail(groupId, "only triangles or tetrahedrons are supported."); + } + + Elements faces = hasTriangles + ? extractSurfaceBoundary(group, coordinates, groupId) + : extractTetrahedronBoundary(group, coordinates, groupId); + const EdgeFaces edgeFaces = buildEdgeFaces(faces, groupId); + validateVertexFans(faces, edgeFaces, groupId); + const auto components = orientComponents(faces, edgeFaces, groupId); + orientOutward(faces, coordinates, components, groupId); + std::sort(faces.begin(), faces.end(), [](const Element& first, const Element& second) { + return faceKey(asFace(first)) < faceKey(asFace(second)); + }); + return faces; +} + +} + +VolumeShellExtractor::VolumeShellExtractor(const Mesh& volumeMesh) : mesh_(volumeMesh) +{ + for (GroupId groupId = 0; groupId < mesh_.groups.size(); ++groupId) { + mesh_.groups[groupId].elements = extractGroupShell( + volumeMesh.groups[groupId], volumeMesh.coordinates, groupId); + } + utils::RedundancyCleaner::cleanCoords(mesh_); +} + +Mesh VolumeShellExtractor::getMesh() const +{ + return mesh_; +} + +} diff --git a/src/core/VolumeShellExtractor.h b/src/core/VolumeShellExtractor.h new file mode 100644 index 0000000..61c85b6 --- /dev/null +++ b/src/core/VolumeShellExtractor.h @@ -0,0 +1,17 @@ +#pragma once + +#include "types/Mesh.h" + +namespace meshlib::core { + +class VolumeShellExtractor { +public: + explicit VolumeShellExtractor(const Mesh& volumeMesh); + + Mesh getMesh() const; + +private: + Mesh mesh_; +}; + +} diff --git a/src/meshers/CMakeLists.txt b/src/meshers/CMakeLists.txt index 32eed77..b48f2ba 100644 --- a/src/meshers/CMakeLists.txt +++ b/src/meshers/CMakeLists.txt @@ -8,12 +8,10 @@ add_library(tessellator-meshers ) target_link_libraries(tessellator-meshers tessellator-core - tessellator-utils - tessellator-cgal - CGAL::CGAL) + tessellator-utils) if(TESSELLATOR_EXECUTION_POLICIES) add_definitions(-DTESSELLATOR_EXECUTION_POLICIES) find_package(TBB CONFIG REQUIRED) target_link_libraries(tessellator-meshers TBB::tbb) -endif() \ No newline at end of file +endif() diff --git a/src/meshers/MesherBase.cpp b/src/meshers/MesherBase.cpp index 20849bb..24c3df0 100644 --- a/src/meshers/MesherBase.cpp +++ b/src/meshers/MesherBase.cpp @@ -52,6 +52,13 @@ void MesherBase::logNumberOfNodes(std::size_t nNodes) log(msg.str(), 2); } +void MesherBase::logNumberOfHexahedra(std::size_t nHexahedra) +{ + std::stringstream msg; + msg << "Mesh contains " << nHexahedra << " hexahedra."; + log(msg.str(), 2); +} + void MesherBase::logGridSize(const Grid& g) { std::stringstream msg; @@ -108,6 +115,7 @@ Mesh MesherBase::buildVolumeMesh(const Mesh& inputMesh, const std::set& Mesh volumeMesh{ inputMesh.grid, inputMesh.coordinates }; volumeMesh.groups.resize(inputMesh.groups.size()); for (const auto& gId : volumeGroups) { + volumeMesh.groups[gId].name = inputMesh.groups[gId].name; mergeGroup(volumeMesh.groups[gId], inputMesh.groups[gId]); } return volumeMesh; @@ -123,4 +131,4 @@ Mesh MesherBase::buildSurfaceMesh(const Mesh& inputMesh, const std::set } } -} \ No newline at end of file +} diff --git a/src/meshers/MesherBase.h b/src/meshers/MesherBase.h index c0bac3e..d05edd2 100644 --- a/src/meshers/MesherBase.h +++ b/src/meshers/MesherBase.h @@ -21,6 +21,7 @@ class MesherBase { static void logNumberOfTriangles(std::size_t nTris); static void logNumberOfLines(std::size_t nLines); static void logNumberOfNodes(std::size_t nNodes); + static void logNumberOfHexahedra(std::size_t nHexahedra); static void logGridSize(const Grid& g); static Grid buildNonSlicingGrid(const Grid& primal, const Grid& enlarged); @@ -37,4 +38,4 @@ class MesherBase { }; } -} \ No newline at end of file +} diff --git a/src/meshers/MesherBaseOptions.h b/src/meshers/MesherBaseOptions.h index 7403619..55138c1 100644 --- a/src/meshers/MesherBaseOptions.h +++ b/src/meshers/MesherBaseOptions.h @@ -7,7 +7,6 @@ namespace meshlib::meshers { class MesherBaseOptions { public: - bool isVolume = false; std::set volumeGroups{}; }; diff --git a/src/meshers/StaircaseMesher.cpp b/src/meshers/StaircaseMesher.cpp index 6e9034d..fadfc34 100644 --- a/src/meshers/StaircaseMesher.cpp +++ b/src/meshers/StaircaseMesher.cpp @@ -7,8 +7,8 @@ #include "core/Collapser.h" #include "core/Staircaser.h" #include "core/Compressor.h" - -#include "cgal/filler/Filler.h" +#include "core/VolumeFiller.h" +#include "core/VolumeShellExtractor.h" #include "utils/RedundancyCleaner.h" #include "utils/MeshTools.h" @@ -20,18 +20,35 @@ using namespace utils; using namespace core; using namespace meshTools; +std::vector getGroupNames(const Groups& groups); +void copyGroupNames(Mesh& mesh, const std::vector& names); + StaircaseMesher::StaircaseMesher(const Mesh& inputMesh, int decimalPlacesInCollapser, StaircaseMesherOptions opts) : MesherBase(inputMesh), decimalPlacesInCollapser_(decimalPlacesInCollapser), opts_(opts) { log("Preparing surfaces."); - surfaceMesh_ = buildMeshFilteringElements(inputMesh, isNotTetrahedron); - + surfaceMesh_ = MesherBase::buildSurfaceMesh(inputMesh, opts_.volumeGroups); log("Processing surface mesh."); process(surfaceMesh_); + + log("Preparing volumes"); + volumeMesh_ = MesherBase::buildVolumeMesh(inputMesh, opts_.volumeGroups); + if (!volumeMesh_.emptyOfElements()) { + volumeMesh_ = VolumeShellExtractor(volumeMesh_).getMesh(); + + log("Processing volume shell."); + process(volumeMesh_, false); + log("Filling volume shell with hexahedra."); + volumeMesh_ = VolumeFiller(volumeMesh_, opts_.splitHexahedra).getMesh(); + logNumberOfHexahedra(countMeshElementsIf(volumeMesh_, isHexahedron)); + } + + mergeMesh(surfaceMesh_, volumeMesh_); + RedundancyCleaner::cleanCoords(surfaceMesh_); - log("Surface mesh built succesfully.", 1); + log("Mesh built succesfully.", 1); } Mesh StaircaseMesher::buildSurfaceMesh(const Mesh& inputMesh, const Mesh & volumeSurface) @@ -41,28 +58,45 @@ Mesh StaircaseMesher::buildSurfaceMesh(const Mesh& inputMesh, const Mesh & volum return resultMesh; } +std::vector getGroupNames(const Groups& groups){ + std::vector names; + names.reserve(groups.size()); + for (auto gId{0}; gId < groups.size(); ++gId) { + names.push_back(groups[gId].name); + } + return names; +} + +void copyGroupNames(Mesh& m, const std::vector& names){ + for (auto gId{0}; gId < m.groups.size(); ++gId) { + m.groups[gId].name = names[gId]; + } +} + +static Mesh toAbsolute(const Mesh& m) +{ + auto r{ m }; + r.coordinates = + utils::GridTools{ m.grid }.relativeToAbsolute(m.coordinates); + return r; +} + void StaircaseMesher::process(Mesh& mesh) const { - + process(mesh, opts_.compress); +} + +void StaircaseMesher::process(Mesh& mesh, bool compress) const +{ + const auto groupNames = getGroupNames(mesh.groups); const auto slicingGrid{ buildSlicingGrid(originalGrid_, enlargedGrid_) }; - if (mesh.countElems() == 0) { - mesh.grid = slicingGrid; + // mesh.grid = slicingGrid; return; } auto dimensions = getHighestDimensionByGroup(mesh); - if (opts_.isVolume){ - if (meshTools::isAClosedTopology(mesh.groups[0].elements)){ - meshlib::cgal::filler::Filler f{ mesh }; - auto filling = f.getMeshFilling(); - mergeMesh(mesh, filling); - } else { - throw std::runtime_error("Input object marked to be meshed as a volume, but surface is not closed"); - } - } - log("Slicing.", 1); mesh.grid = slicingGrid; mesh = Slicer{ mesh, dimensions }.getMesh(); @@ -86,7 +120,7 @@ void StaircaseMesher::process(Mesh& mesh) const logNumberOfQuads(countMeshElementsIf(mesh, isQuad)); logNumberOfLines(countMeshElementsIf(mesh, isLine)); - if (opts_.compress) { + if (compress) { log("Compressing surfaces.", 1); std::size_t beforeQuads = countMeshElementsIf(mesh, isQuad); std::size_t merged = Compressor::compressSurfacesInMesh(mesh); @@ -113,6 +147,8 @@ void StaircaseMesher::process(Mesh& mesh) const logNumberOfQuads(countMeshElementsIf(mesh, isQuad)); logNumberOfLines(countMeshElementsIf(mesh, isLine)); + copyGroupNames(mesh, groupNames); + } @@ -121,4 +157,4 @@ Mesh StaircaseMesher::mesh() const return surfaceMesh_; } -} \ No newline at end of file +} diff --git a/src/meshers/StaircaseMesher.h b/src/meshers/StaircaseMesher.h index 3cabd4b..d768399 100644 --- a/src/meshers/StaircaseMesher.h +++ b/src/meshers/StaircaseMesher.h @@ -17,10 +17,12 @@ class StaircaseMesher : public MesherBase { int decimalPlacesInCollapser_; Mesh surfaceMesh_; + Mesh volumeMesh_; StaircaseMesherOptions opts_; virtual Mesh buildSurfaceMesh(const Mesh& inputMesh, const Mesh& volumeSurface); void process(Mesh&) const; + void process(Mesh&, bool compress) const; }; diff --git a/src/meshers/StaircaseMesherOptions.h b/src/meshers/StaircaseMesherOptions.h index 53f2b50..4bac022 100644 --- a/src/meshers/StaircaseMesherOptions.h +++ b/src/meshers/StaircaseMesherOptions.h @@ -7,6 +7,7 @@ namespace meshlib::meshers { class StaircaseMesherOptions : public MesherBaseOptions { public: bool compress = false; + bool splitHexahedra = false; }; } diff --git a/src/types/Mesh.h b/src/types/Mesh.h index dcccb37..c96d3b8 100644 --- a/src/types/Mesh.h +++ b/src/types/Mesh.h @@ -69,6 +69,11 @@ struct Element { return type == Type::Volume && vertices.size() == 4; } + bool isHexahedron() const + { + return type == Type::Volume && vertices.size() == 8; + } + bool sharesVertices(const Element& rhs) { bool res = true; @@ -220,4 +225,3 @@ struct Mesh { }; } - diff --git a/src/utils/Geometry.cpp b/src/utils/Geometry.cpp index fe2fc4e..2015763 100644 --- a/src/utils/Geometry.cpp +++ b/src/utils/Geometry.cpp @@ -70,6 +70,17 @@ std::vector Geometry::buildDisjointSmoothSets( } +QuaV Geometry::asQuaV(const Element& el, const std::vector& co) { + if (el.vertices.size() != 4) { + throw std::logic_error("Invalid conversion from element to QuaV"); + } + QuaV res; + for (std::size_t i = 0; i < el.vertices.size(); i++) { + res[i] = co[el.vertices[i]]; + } + return res; +} + TriV Geometry::asTriV(const Element& el, const std::vector& co) { if (el.vertices.size() != 3) { throw std::logic_error("Invalid conversion from element to TriV"); @@ -229,6 +240,18 @@ double Geometry::area(const TriV& tri) { return ((tri[0] - tri[1]) ^ (tri[1] - tri[2])).norm() / 2.0; } +double Geometry::area(const QuaV& qua) { + const Coordinates cs{ qua.begin(), qua.end() }; + + VecD crossSum; + for (std::size_t i = 0; i < qua.size(); ++i) { + const auto& p = qua[i]; + const auto& q = qua[(i + 1) % qua.size()]; + crossSum += (p ^ q); + } + return 0.5*crossSum.norm(); +} + } } diff --git a/src/utils/Geometry.h b/src/utils/Geometry.h index 113593b..e185f0c 100644 --- a/src/utils/Geometry.h +++ b/src/utils/Geometry.h @@ -28,6 +28,7 @@ class Geometry { static bool areAdjacentLines(const Element&, const Element&); + static QuaV asQuaV(const Element&, const Coordinates&); static TriV asTriV(const Element&, const Coordinates&); static LinV asLinV(const Element&, const Coordinates&); @@ -38,6 +39,7 @@ class Geometry { static VecD getCentroid(const Element&, const std::vector&); static VecD getCentroid(const TriV&); static double area(const TriV& tri); + static double area(const QuaV& qua); static bool isDegenerate(const TriV& tri, const double& areaTolerance = NORM_TOLERANCE); static bool areCollinear(const Coordinates&); template diff --git a/src/utils/MeshTools.cpp b/src/utils/MeshTools.cpp index 02bbe57..00e0f1b 100644 --- a/src/utils/MeshTools.cpp +++ b/src/utils/MeshTools.cpp @@ -273,12 +273,22 @@ void checkNoNullAreasExist(const Mesh& m) msg << info(e, m) << std::endl; } } - else if (Geometry::area(Geometry::asTriV(e, m.coordinates)) == 0.0) { - nullAreas = true; - msg << std::endl; - msg << "Group: " << &g - &m.groups.front() - << ", Element: " << &e - &g.elements.front() << std::endl; - msg << info(e, m) << std::endl; + else if (e.isTriangle()){ + if (Geometry::area(Geometry::asTriV(e, m.coordinates)) == 0.0) { + nullAreas = true; + msg << std::endl; + msg << "Group: " << &g - &m.groups.front() + << ", Element: " << &e - &g.elements.front() << std::endl; + msg << info(e, m) << std::endl; + } + } else if (e.isQuad()){ + if (Geometry::area(Geometry::asQuaV(e, m.coordinates)) == 0.0) { + nullAreas = true; + msg << std::endl; + msg << "Group: " << &g - &m.groups.front() + << ", Element: " << &e - &g.elements.front() << std::endl; + msg << info(e, m) << std::endl; + } } } } @@ -299,6 +309,17 @@ void convertToAbsoluteCoordinates(Mesh& m) ); } +void convertToRelativeCoordinates(Mesh& m) +{ + GridTools gT{ m.grid }; + + std::transform( + m.coordinates.begin(), m.coordinates.end(), + m.coordinates.begin(), + [&](const auto& v) { return gT.getRelative(v); } + ); +} + void checkSlicedMeshInvariants(const Mesh& m) { checkNoCellsAreCrossed(m); @@ -319,6 +340,7 @@ Mesh buildMeshFilteringElements( inElems.begin(), inElems.end(), std::back_inserter(r.groups[gId].elements), filter); + r.groups[gId].name = in.groups[gId].name; } return r; } @@ -378,7 +400,7 @@ void mergeMesh(Mesh& lMesh, const Mesh& iMesh) assert(lMesh.groups.size() == iMesh.groups.size()); auto coordCount{ lMesh.coordinates.size() }; - + if (iMesh.countElems() == 0) return; lMesh.coordinates.insert(lMesh.coordinates.end(), iMesh.coordinates.begin(), iMesh.coordinates.end()); diff --git a/src/utils/MeshTools.h b/src/utils/MeshTools.h index 5fc4546..4826a0d 100644 --- a/src/utils/MeshTools.h +++ b/src/utils/MeshTools.h @@ -19,6 +19,8 @@ static bool isQuad(const Element& e) { return e.isQuad(); } static bool isNotQuad(const Element& e) { return !e.isQuad(); } static bool isTetrahedron(const Element& e) { return e.isTetrahedron(); } static bool isNotTetrahedron(const Element& e) { return !e.isTetrahedron(); } +static bool isHexahedron(const Element& e) { return e.isHexahedron(); } +static bool isNotHexahedron(const Element& e) { return !e.isHexahedron(); } std::size_t countMeshElementsIf(const Mesh& mesh, std::function countFilter); std::vector getHighestDimensionByGroup(const Mesh& mesh); @@ -35,6 +37,7 @@ void reduceGrid(Mesh&, const Grid&); Mesh reduceGrid(const Mesh& m, const Grid& g); void convertToAbsoluteCoordinates(Mesh&); +void convertToRelativeCoordinates(Mesh&); void checkSlicedMeshInvariants(Mesh& m); @@ -52,4 +55,4 @@ bool isAClosedTopology(const Elements& es); Mesh extractGroupsByName(const Mesh& mesh, const std::vector& groupNames); -} \ No newline at end of file +} diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 638a03a..3cc1202 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -20,7 +20,9 @@ add_executable(tessellator_tests "core/SmootherTest.cpp" "core/SmootherToolsTest.cpp" "core/SnapperTest.cpp" - "core/StaircaserTest.cpp" + "core/StaircaserTest.cpp" + "core/VolumeFillerTest.cpp" + "core/VolumeShellExtractorTest.cpp" "types/MeshTest.cpp" "utils/ConvexHullTest.cpp" "utils/CoordGraphTest.cpp" @@ -67,4 +69,4 @@ if (TESSELLATOR_ENABLE_CGAL) target_link_libraries(tessellator_tests tessellator-cgal) endif() -add_test(tessellator ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/tessellator_tests) \ No newline at end of file +add_test(tessellator ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/tessellator_tests) diff --git a/test/MeshFixtures.h b/test/MeshFixtures.h index 65c77ad..7eda7c0 100644 --- a/test/MeshFixtures.h +++ b/test/MeshFixtures.h @@ -38,6 +38,7 @@ static Mesh buildNonManifoldPatchMesh(double stepSize) return m; } + static Mesh buildTetAndTriMesh(double stepSize) { Mesh m; diff --git a/test/app/launcherTest.cpp b/test/app/launcherTest.cpp index 5b3f3e4..2cf01ab 100644 --- a/test/app/launcherTest.cpp +++ b/test/app/launcherTest.cpp @@ -71,8 +71,32 @@ TEST_F(LauncherTest, builds_staircased_mesher_default) EXPECT_NO_THROW(auto staircaseMesher = dynamic_cast(*mesher)); const auto & options = dynamic_cast(*mesher).getOptions(); - EXPECT_EQ(options.isVolume, false); + EXPECT_EQ(options.volumeGroups.size(), 0); EXPECT_EQ(options.compress, false); + EXPECT_FALSE(options.splitHexahedra); +} + +TEST_F(LauncherTest, buildsStaircasedMesherWithSplitHexahedra) +{ + meshlib::Mesh meshMock; + meshMock.grid = { + std::vector{0, 1}, + std::vector{0, 1}, + std::vector{0, 1} + }; + nlohmann::json config = { + {"mesher", { + {"type", "staircase"}, + {"options", {{"splitHexahedra", true}}} + }} + }; + + ObjectDefinition object; + auto mesher = buildMesher(meshMock, config, object); + const auto& staircase = + dynamic_cast(*mesher); + + EXPECT_TRUE(staircase.getOptions().splitHexahedra); } TEST_F(LauncherTest, builds_staircased_mesher_without_compression) @@ -95,7 +119,7 @@ TEST_F(LauncherTest, builds_staircased_mesher_without_compression) EXPECT_NO_THROW(auto staircaseMesher = dynamic_cast(*mesher)); const auto & options = dynamic_cast(*mesher).getOptions(); - EXPECT_EQ(options.isVolume, false); + EXPECT_EQ(options.volumeGroups.size(), 0); EXPECT_EQ(options.compress, false); } @@ -119,7 +143,7 @@ TEST_F(LauncherTest, builds_staircased_mesher_with_compression) EXPECT_NO_THROW(auto staircaseMesher = dynamic_cast(*mesher)); const auto & options = dynamic_cast(*mesher).getOptions(); - EXPECT_EQ(options.isVolume, false); + EXPECT_EQ(options.volumeGroups.size(), 0); EXPECT_EQ(options.compress, true); } @@ -315,7 +339,7 @@ TEST_F(LauncherTest, builds_staircased_mesher_with_override) EXPECT_NO_THROW(auto staircaseMesher = dynamic_cast(*mesher)); const auto & options = dynamic_cast(*mesher).getOptions(); - EXPECT_EQ(options.isVolume, false); + EXPECT_EQ(options.volumeGroups.size(), 0); EXPECT_EQ(options.compress, true); } @@ -354,4 +378,3 @@ TEST_F(LauncherTest, launches_multiObject_sameFileMultipleGroups) 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 5541b85..4e89dc5 100644 --- a/test/app/vtkIOTest.cpp +++ b/test/app/vtkIOTest.cpp @@ -43,6 +43,30 @@ TEST_F(VTKIOTest, readElementTypes) EXPECT_TRUE(m.groups[0].elements[2].isTriangle()); } +TEST_F(VTKIOTest, exportAndReadHexahedron) +{ + meshlib::Mesh mesh; + mesh.grid = meshlib::utils::GridTools::buildCartesianGrid(0.0, 1.0, 2); + mesh.coordinates = { + meshlib::Coordinate({0.0, 0.0, 0.0}), meshlib::Coordinate({1.0, 0.0, 0.0}), + meshlib::Coordinate({1.0, 1.0, 0.0}), meshlib::Coordinate({0.0, 1.0, 0.0}), + meshlib::Coordinate({0.0, 0.0, 1.0}), meshlib::Coordinate({1.0, 0.0, 1.0}), + meshlib::Coordinate({1.0, 1.0, 1.0}), meshlib::Coordinate({0.0, 1.0, 1.0})}; + mesh.groups.resize(1); + mesh.groups[0].elements.push_back(meshlib::Element( + {0, 1, 2, 3, 4, 5, 6, 7}, meshlib::Element::Type::Volume)); + + const auto filename = std::filesystem::temp_directory_path() + / "tessellator_hexahedron_roundtrip.vtu"; + exportMeshToVTU(filename, mesh); + const auto result = readInputMesh(filename); + std::filesystem::remove(filename); + + ASSERT_EQ(1, result.countElems()); + EXPECT_TRUE(result.groups[0].elements[0].isHexahedron()); + EXPECT_EQ(mesh.coordinates, result.coordinates); +} + TEST_F(VTKIOTest, exportGridToVTU) { meshlib::Grid grid; @@ -56,4 +80,4 @@ TEST_F(VTKIOTest, exportGridToVTU) auto exported{ readInputMesh(fn) }; EXPECT_EQ(121+121+21, exported.countElems()); -} \ No newline at end of file +} diff --git a/test/cgal/ManifolderTest.cpp b/test/cgal/ManifolderTest.cpp index cedd1c3..c05f1c5 100644 --- a/test/cgal/ManifolderTest.cpp +++ b/test/cgal/ManifolderTest.cpp @@ -87,6 +87,7 @@ TEST_F(ManifolderTest, volume_and_surface) ASSERT_EQ(1, r.countElems()); } + TEST_F(ManifolderTest, closed_surface) { Mesh m = buildCubeSurfaceMesh(1.0); diff --git a/test/cgal/filler/FillerTest.cpp b/test/cgal/filler/FillerTest.cpp index b2230e6..e8fbfa7 100644 --- a/test/cgal/filler/FillerTest.cpp +++ b/test/cgal/filler/FillerTest.cpp @@ -155,6 +155,14 @@ class FillerTest : public ::testing::Test { return r; } + static Mesh toAbsolute(const Mesh& m) + { + auto r{ m }; + r.coordinates = + utils::GridTools{ m.grid }.relativeToAbsolute(m.coordinates); + return r; + } + static bool allAreSimple(const FaceFilling& ff) { for (const auto [pr, ss] : ff.tris) { @@ -337,6 +345,59 @@ TEST_F(FillerTest, parallelogram_as_surface) EXPECT_EQ(0, countPWHs(f.getFaceFilling({ Cell({0, 0, 1}), Z }))); } +TEST_F(FillerTest, fill_cube1x1x1_size1_grid) +{ + Mesh m = buildCubeSurfaceMesh(01.0); + + Mesh out; + ASSERT_NO_THROW(out = Slicer{m}.getMesh()); + EXPECT_EQ(12, countMeshElementsIf(out, isTriangle)); + + Mesh filled = Filler{out}.getMeshFilling(); + EXPECT_EQ(12, countMeshElementsIf(filled, isTriangle)); +} + +TEST_F(FillerTest, fill_cube1x1x1_size05_grid) +{ + //filling with unstruc. triangles + Mesh m = buildCubeSurfaceMesh(0.5); + Mesh filled = Filler{Slicer{buildCubeSurfaceMesh(0.5) }.getMesh()}.getMeshFilling(); + Mesh filled_no_slicing = Filler{toRelative(buildCubeSurfaceMesh(0.5))}.getMeshFilling(); + EXPECT_EQ(18, countMeshElementsIf(filled, isTriangle)); + EXPECT_EQ(18, countMeshElementsIf(filled_no_slicing, isTriangle)); + + //slicing hull + Mesh out_1; + ASSERT_NO_THROW(out_1 = Slicer{buildCubeSurfaceMesh(0.5) }.getMesh()); + EXPECT_EQ(48, countMeshElementsIf(out_1, isTriangle)); + + //filling and slicing + Mesh out_2; + ASSERT_NO_THROW(out_2 = Slicer{toAbsolute(filled) }.getMesh()); + EXPECT_EQ(72, countMeshElementsIf(out_2, isTriangle)); + + ASSERT_NO_THROW(out_2 = Slicer{toAbsolute(filled_no_slicing)}.getMesh()); + EXPECT_EQ(72, countMeshElementsIf(out_2, isTriangle)); + +} + +TEST_F(FillerTest, fill_cube1x1x1_size025_grid) +{ + Mesh m = buildCubeSurfaceMesh(0.25); + Mesh filled = Filler{Slicer{buildCubeSurfaceMesh(0.25) }.getMesh()}.getMeshFilling(); + EXPECT_EQ(30, countMeshElementsIf(filled, isTriangle)); + + Mesh out_1; + ASSERT_NO_THROW(out_1 = Slicer{buildCubeSurfaceMesh(0.25) }.getMesh()); + EXPECT_EQ(192, countMeshElementsIf(out_1, isTriangle)); + + Mesh out_2; + ASSERT_NO_THROW(out_2 = Slicer{toAbsolute(filled) }.getMesh()); + EXPECT_EQ(480, countMeshElementsIf(out_2, isTriangle)); + +} + + TEST_F(FillerTest, planeXY_mesh_filling) { Filler f{ Slicer{ buildPlaneXYMesh(1.0) }.getMesh() }; diff --git a/test/core/VolumeFillerTest.cpp b/test/core/VolumeFillerTest.cpp new file mode 100644 index 0000000..01b7fa8 --- /dev/null +++ b/test/core/VolumeFillerTest.cpp @@ -0,0 +1,148 @@ +#include "gtest/gtest.h" + +#include "core/VolumeFiller.h" +#include "utils/GridTools.h" +#include "utils/MeshTools.h" + +namespace meshlib::core { + +using namespace utils; +using namespace meshTools; + +namespace { + +CoordinateId coordinateId(Mesh& mesh, CellDir x, CellDir y, CellDir z) +{ + const Coordinate coordinate = GridTools(mesh.grid).getPos(Cell({x, y, z})); + const auto found = std::find(mesh.coordinates.begin(), mesh.coordinates.end(), coordinate); + if (found != mesh.coordinates.end()) { + return found - mesh.coordinates.begin(); + } + mesh.coordinates.push_back(coordinate); + return mesh.coordinates.size() - 1; +} + +void addQuad(Mesh& mesh, const std::array& cells) +{ + Element quad; + quad.type = Element::Type::Surface; + for (const auto& cell : cells) { + quad.vertices.push_back(coordinateId(mesh, cell[X], cell[Y], cell[Z])); + } + mesh.groups[0].elements.push_back(quad); +} + +Mesh buildTwoByTwoByTwoShell() +{ + Mesh mesh; + mesh.grid = GridTools::buildCartesianGrid(0.0, 1.0, 3); + mesh.groups.resize(1); + mesh.groups[0].name = "volume"; + + for (CellDir first = 0; first < 2; ++first) { + for (CellDir second = 0; second < 2; ++second) { + for (CellDir x : {CellDir(0), CellDir(2)}) { + addQuad(mesh, { + Cell({x, first, second}), Cell({x, first + 1, second}), + Cell({x, first + 1, second + 1}), Cell({x, first, second + 1})}); + } + for (CellDir y : {CellDir(0), CellDir(2)}) { + addQuad(mesh, { + Cell({first, y, second}), Cell({first + 1, y, second}), + Cell({first + 1, y, second + 1}), Cell({first, y, second + 1})}); + } + for (CellDir z : {CellDir(0), CellDir(2)}) { + addQuad(mesh, { + Cell({first, second, z}), Cell({first + 1, second, z}), + Cell({first + 1, second + 1, z}), Cell({first, second + 1, z})}); + } + } + } + return mesh; +} + +void addBox(Mesh& mesh, CellDir x0, CellDir x1) +{ + addQuad(mesh, {Cell({x0, 0, 0}), Cell({x0, 1, 0}), + Cell({x0, 1, 1}), Cell({x0, 0, 1})}); + addQuad(mesh, {Cell({x1, 0, 0}), Cell({x1, 1, 0}), + Cell({x1, 1, 1}), Cell({x1, 0, 1})}); + addQuad(mesh, {Cell({x0, 0, 0}), Cell({x1, 0, 0}), + Cell({x1, 0, 1}), Cell({x0, 0, 1})}); + addQuad(mesh, {Cell({x0, 1, 0}), Cell({x1, 1, 0}), + Cell({x1, 1, 1}), Cell({x0, 1, 1})}); + addQuad(mesh, {Cell({x0, 0, 0}), Cell({x1, 0, 0}), + Cell({x1, 1, 0}), Cell({x0, 1, 0})}); + addQuad(mesh, {Cell({x0, 0, 1}), Cell({x1, 0, 1}), + Cell({x1, 1, 1}), Cell({x0, 1, 1})}); +} + +} + +TEST(VolumeFillerTest, fillsContinuousRunsWithHexahedra) +{ + const Mesh result = VolumeFiller(buildTwoByTwoByTwoShell()).getMesh(); + + EXPECT_EQ("volume", result.groups[0].name); + EXPECT_EQ(4, countMeshElementsIf(result, isHexahedron)); + EXPECT_EQ(0, countMeshElementsIf(result, isQuad)); + EXPECT_EQ(18, result.coordinates.size()); + + const auto& first = result.groups[0].elements.front(); + const Coordinates expected = { + Coordinate({0.0, 0.0, 0.0}), Coordinate({1.0, 0.0, 0.0}), + Coordinate({1.0, 0.5, 0.0}), Coordinate({0.0, 0.5, 0.0}), + Coordinate({0.0, 0.0, 0.5}), Coordinate({1.0, 0.0, 0.5}), + Coordinate({1.0, 0.5, 0.5}), Coordinate({0.0, 0.5, 0.5})}; + for (std::size_t vertex = 0; vertex < expected.size(); ++vertex) { + EXPECT_EQ(expected[vertex], result.coordinates[first.vertices[vertex]]); + } +} + +TEST(VolumeFillerTest, splitsContinuousRunsIntoUnitCellHexahedra) +{ + const Mesh result = VolumeFiller(buildTwoByTwoByTwoShell(), true).getMesh(); + + EXPECT_EQ(8, countMeshElementsIf(result, isHexahedron)); + const GridTools tools(result.grid); + for (const auto& element : result.groups[0].elements) { + ASSERT_TRUE(element.isHexahedron()); + Cell lower = tools.getCell(result.coordinates[element.vertices[0]]); + Cell upper = lower; + for (CoordinateId vertex : element.vertices) { + const Cell cell = tools.getCell(result.coordinates[vertex]); + for (Axis axis : {X, Y, Z}) { + lower[axis] = std::min(lower[axis], cell[axis]); + upper[axis] = std::max(upper[axis], cell[axis]); + } + } + for (Axis axis : {X, Y, Z}) { + EXPECT_EQ(1, upper[axis] - lower[axis]); + } + } +} + +TEST(VolumeFillerTest, rejectsAnOpenQuadShell) +{ + Mesh shell = buildTwoByTwoByTwoShell(); + shell.groups[0].elements.erase(shell.groups[0].elements.begin()); + + EXPECT_THROW(VolumeFiller{shell}, std::runtime_error); +} + +TEST(VolumeFillerTest, fillsDisconnectedIntervalsOnTheSameRay) +{ + Mesh shell; + shell.grid = GridTools::buildCartesianGrid(0.0, 3.0, 4); + shell.groups.resize(1); + addBox(shell, 0, 1); + addBox(shell, 2, 3); + + const Mesh result = VolumeFiller(shell).getMesh(); + + ASSERT_EQ(2, result.countElems()); + EXPECT_TRUE(result.groups[0].elements[0].isHexahedron()); + EXPECT_TRUE(result.groups[0].elements[1].isHexahedron()); +} + +} diff --git a/test/core/VolumeShellExtractorTest.cpp b/test/core/VolumeShellExtractorTest.cpp new file mode 100644 index 0000000..8b9d656 --- /dev/null +++ b/test/core/VolumeShellExtractorTest.cpp @@ -0,0 +1,175 @@ +#include "gtest/gtest.h" + +#include "MeshFixtures.h" +#include "core/VolumeShellExtractor.h" +#include "utils/Geometry.h" +#include "utils/MeshTools.h" + +namespace meshlib::core { + +using namespace meshFixtures; +using namespace utils::meshTools; + +namespace { + +double signedVolume6(const Mesh& mesh, const Group& group) +{ + double result = 0.0; + for (const Element& element : group.elements) { + const auto& first = mesh.coordinates[element.vertices[0]]; + const auto& second = mesh.coordinates[element.vertices[1]]; + const auto& third = mesh.coordinates[element.vertices[2]]; + result += first * (second ^ third); + } + return result; +} + +Mesh buildOpenTetrahedronShell() +{ + Mesh mesh = buildTetSurfaceMesh(1.0); + mesh.groups[0].elements.pop_back(); + return mesh; +} + +} + +TEST(VolumeShellExtractorTest, extractsOutwardBoundaryFromTetrahedron) +{ + const Mesh shell = VolumeShellExtractor(buildTetMesh(1.0)).getMesh(); + + EXPECT_EQ(4, countMeshElementsIf(shell, isTriangle)); + EXPECT_TRUE(isAClosedTopology(shell.groups[0].elements)); + EXPECT_GT(signedVolume6(shell, shell.groups[0]), 0.0); +} + +TEST(VolumeShellExtractorTest, removesFacesInsideConnectedTetrahedrons) +{ + const Mesh shell = VolumeShellExtractor(buildTetMeshWithInnerPoint(1.0)).getMesh(); + + EXPECT_EQ(4, countMeshElementsIf(shell, isTriangle)); + EXPECT_EQ(4, shell.coordinates.size()); + EXPECT_TRUE(isAClosedTopology(shell.groups[0].elements)); +} + +TEST(VolumeShellExtractorTest, acceptsAndOrientsClosedTriangleShell) +{ + Mesh input = buildTetSurfaceMesh(1.0); + std::reverse( + input.groups[0].elements[0].vertices.begin(), + input.groups[0].elements[0].vertices.end()); + + const Mesh shell = VolumeShellExtractor(input).getMesh(); + + EXPECT_EQ(4, shell.countElems()); + EXPECT_TRUE(isAClosedTopology(shell.groups[0].elements)); + EXPECT_GT(signedVolume6(shell, shell.groups[0]), 0.0); +} + +TEST(VolumeShellExtractorTest, preservesEmptyGroupsAndNames) +{ + Mesh input = buildTetMesh(1.0); + input.groups[0].name = "volume"; + input.groups.insert(input.groups.begin(), Group{"surface", {}}); + + const Mesh shell = VolumeShellExtractor(input).getMesh(); + + ASSERT_EQ(2, shell.groups.size()); + EXPECT_EQ("surface", shell.groups[0].name); + EXPECT_TRUE(shell.groups[0].elements.empty()); + EXPECT_EQ("volume", shell.groups[1].name); + EXPECT_EQ(4, shell.groups[1].elements.size()); +} + +TEST(VolumeShellExtractorTest, acceptsDisconnectedClosedComponents) +{ + Mesh input = buildTetMesh(1.0); + input.coordinates.insert(input.coordinates.end(), { + Coordinate({2.0, 0.0, 0.0}), + Coordinate({3.0, 0.0, 0.0}), + Coordinate({2.0, 1.0, 0.0}), + Coordinate({2.0, 0.0, 1.0})}); + input.groups[0].elements.push_back( + Element({4, 5, 6, 7}, Element::Type::Volume)); + + const Mesh shell = VolumeShellExtractor(input).getMesh(); + + EXPECT_EQ(8, shell.countElems()); + EXPECT_GT(signedVolume6(shell, shell.groups[0]), 0.0); +} + +TEST(VolumeShellExtractorTest, rejectsOpenShell) +{ + EXPECT_THROW(VolumeShellExtractor{buildOpenTetrahedronShell()}, std::runtime_error); +} + +TEST(VolumeShellExtractorTest, rejectsMixedTrianglesAndTetrahedrons) +{ + Mesh input = buildTetMesh(1.0); + input.groups[0].elements.push_back(Element({0, 1, 2})); + + EXPECT_THROW(VolumeShellExtractor{input}, std::runtime_error); +} + +TEST(VolumeShellExtractorTest, rejectsNonManifoldEdge) +{ + Mesh input = buildTetSurfaceMesh(1.0); + input.coordinates.push_back(Coordinate({0.0, -1.0, 0.0})); + input.groups[0].elements.push_back(Element({0, 1, 4})); + + EXPECT_THROW(VolumeShellExtractor{input}, std::runtime_error); +} + +TEST(VolumeShellExtractorTest, rejectsTetrahedronsTouchingOnlyAtOneVertex) +{ + Mesh input = buildTetMesh(1.0); + input.coordinates.insert(input.coordinates.end(), { + Coordinate({-1.0, 0.0, 0.0}), + Coordinate({0.0, -1.0, 0.0}), + Coordinate({0.0, 0.0, -1.0})}); + input.groups[0].elements.push_back( + Element({0, 4, 5, 6}, Element::Type::Volume)); + + EXPECT_THROW(VolumeShellExtractor{input}, std::runtime_error); +} + +TEST(VolumeShellExtractorTest, rejectsTetrahedronsTouchingOnlyAtOneEdge) +{ + EXPECT_THROW( + VolumeShellExtractor{buildTetsSharingEdgeMesh()}, + std::runtime_error); +} + +TEST(VolumeShellExtractorTest, rejectsFaceSharedByThreeTetrahedrons) +{ + Mesh input = buildTetMesh(1.0); + input.coordinates.push_back(Coordinate({0.0, 0.0, -1.0})); + input.coordinates.push_back(Coordinate({0.2, 0.2, -1.0})); + input.groups[0].elements.push_back( + Element({0, 2, 1, 4}, Element::Type::Volume)); + input.groups[0].elements.push_back( + Element({0, 1, 2, 5}, Element::Type::Volume)); + + EXPECT_THROW(VolumeShellExtractor{input}, std::runtime_error); +} + +TEST(VolumeShellExtractorTest, rejectsInvalidCoordinateId) +{ + Mesh input = buildTetMesh(1.0); + input.groups[0].elements[0].vertices[3] = input.coordinates.size(); + + EXPECT_THROW(VolumeShellExtractor{input}, std::runtime_error); +} + +TEST(VolumeShellExtractorTest, rejectsDegenerateAndUnsupportedElements) +{ + Mesh degenerate = buildTetMesh(1.0); + degenerate.coordinates[3] = Coordinate({0.5, 0.5, 0.0}); + EXPECT_THROW(VolumeShellExtractor{degenerate}, std::runtime_error); + + Mesh unsupported = buildTetMesh(1.0); + unsupported.groups[0].elements = { + Element({0, 1, 2, 3, 0, 1, 2, 3}, Element::Type::Volume)}; + EXPECT_THROW(VolumeShellExtractor{unsupported}, std::runtime_error); +} + +} diff --git a/test/meshers/StaircaseMesherTest.cpp b/test/meshers/StaircaseMesherTest.cpp index 6d1dfdb..2ec23d1 100644 --- a/test/meshers/StaircaseMesherTest.cpp +++ b/test/meshers/StaircaseMesherTest.cpp @@ -1,6 +1,10 @@ #include "gtest/gtest.h" #include "MeshFixtures.h" +#include +#include +#include +#include #include "meshers/StaircaseMesher.h" #include "StaircaseMesherOptions.h" @@ -24,6 +28,132 @@ using namespace meshFixtures; using namespace utils; using namespace meshTools; +namespace { + +using EdgeKey = std::array; +using FaceKey = std::array; + +EdgeKey edgeKey(CoordinateId first, CoordinateId second) +{ + return first < second ? EdgeKey{first, second} : EdgeKey{second, first}; +} + +bool isSingleClosedSurface(const Elements& faces) +{ + if (faces.empty()) { + return false; + } + std::map> edgeFaces; + for (ElementId faceId = 0; faceId < faces.size(); ++faceId) { + const auto& vertices = faces[faceId].vertices; + if (vertices.size() < 3) { + return false; + } + for (std::size_t vertex = 0; vertex < vertices.size(); ++vertex) { + edgeFaces[edgeKey(vertices[vertex], vertices[(vertex + 1) % vertices.size()])] + .push_back(faceId); + } + } + + std::vector> adjacency(faces.size()); + for (const auto& entry : edgeFaces) { + if (entry.second.size() != 2) { + return false; + } + const ElementId first = entry.second[0]; + const ElementId second = entry.second[1]; + adjacency[first].insert(second); + adjacency[second].insert(first); + } + + std::set visited{0}; + std::queue pending; + pending.push(0); + while (!pending.empty()) { + const ElementId current = pending.front(); + pending.pop(); + for (ElementId neighbor : adjacency[current]) { + if (visited.insert(neighbor).second) { + pending.push(neighbor); + } + } + } + return visited.size() == faces.size(); +} + +bool isSingleClosedHexahedralVolume(const Mesh& mesh) +{ + static const std::array, 6> hexahedronFaces{{ + {{0, 1, 2, 3}}, {{4, 5, 6, 7}}, + {{0, 1, 5, 4}}, {{1, 2, 6, 5}}, + {{2, 3, 7, 6}}, {{3, 0, 4, 7}} + }}; + + std::vector hexahedra; + for (const auto& group : mesh.groups) { + for (const auto& element : group.elements) { + if (!element.isHexahedron()) { + return false; + } + hexahedra.push_back(&element); + } + } + if (hexahedra.empty()) { + return false; + } + + struct FaceOccurrence { + ElementId hexahedron; + CoordinateIds vertices; + }; + std::map> faceOccurrences; + for (ElementId elementId = 0; elementId < hexahedra.size(); ++elementId) { + for (const auto& face : hexahedronFaces) { + CoordinateIds vertices; + FaceKey key; + for (std::size_t vertex = 0; vertex < face.size(); ++vertex) { + vertices.push_back(hexahedra[elementId]->vertices[face[vertex]]); + key[vertex] = vertices.back(); + } + std::sort(key.begin(), key.end()); + faceOccurrences[key].push_back({elementId, vertices}); + } + } + + Elements boundary; + std::vector> adjacency(hexahedra.size()); + for (const auto& entry : faceOccurrences) { + if (entry.second.size() == 1) { + boundary.emplace_back( + entry.second.front().vertices, Element::Type::Surface); + } else if (entry.second.size() == 2) { + const ElementId first = entry.second[0].hexahedron; + const ElementId second = entry.second[1].hexahedron; + adjacency[first].insert(second); + adjacency[second].insert(first); + } else { + return false; + } + } + + std::set visited{0}; + std::queue pending; + pending.push(0); + while (!pending.empty()) { + const ElementId current = pending.front(); + pending.pop(); + for (ElementId neighbor : adjacency[current]) { + if (visited.insert(neighbor).second) { + pending.push(neighbor); + } + } + } + return visited.size() == hexahedra.size() + && isSingleClosedSurface(boundary); +} + +} + class StaircaseMesherTest : public ::testing::Test { public: @@ -369,9 +499,84 @@ TEST_F(StaircaseMesherTest, testStaircaseWithCompression) EXPECT_EQ(0, countMeshElementsIf(compressedMesh, isNode)); } +TEST_F(StaircaseMesherTest, mesh_tetrahedron_volume_2x2){ + + Mesh m = buildCubeVolumeMesh(0.5); + meshlib::meshers::StaircaseMesherOptions opts; + opts.volumeGroups.insert(0); + +// #if APP_LOADED +// vtkIO::exportMeshToVTU("testData/cases/mesh_tetrahedron_volume_2x2_before.vtk", m); +// vtkIO::exportGridToVTU("testData/cases/mesh_tetrahedron_volume_2x2_before_grid.vtk", m.grid); +// #endif + + auto staircasedMesh = StaircaseMesher{m, 4, opts }.mesh(); + +// #if APP_LOADED +// vtkIO::exportMeshToVTU("testData/cases/mesh_tetrahedron_volume_2x2_after.vtk", staircasedMesh); +// vtkIO::exportGridToVTU("testData/cases/mesh_tetrahedron_volume_2x2_after_grid.vtk", staircasedMesh.grid); +// #endif + + EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isTriangle)); + EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isQuad)); + EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isTetrahedron)); + EXPECT_EQ(4, countMeshElementsIf(staircasedMesh, isHexahedron)); + EXPECT_EQ(18, staircasedMesh.coordinates.size()); + +} + +TEST_F(StaircaseMesherTest, mesh_surface_volume_2x2){ + + Mesh m = buildCubeSurfaceMesh(0.5); + meshlib::meshers::StaircaseMesherOptions opts; + opts.volumeGroups.insert(0); + // opts.isVolume = true; + auto staircasedMesh = StaircaseMesher{m, 4, opts }.mesh(); + + EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isTriangle)); + EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isQuad)); + EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isTetrahedron)); + EXPECT_EQ(4, countMeshElementsIf(staircasedMesh, isHexahedron)); + EXPECT_EQ(18, staircasedMesh.coordinates.size()); + +} + +TEST_F(StaircaseMesherTest, mesh_surface_not_volume_2x2){ + + Mesh m = buildCubeSurfaceMesh(0.5); + meshlib::meshers::StaircaseMesherOptions opts; + // opts.isVolume = true; + auto staircasedMesh = StaircaseMesher{m, 4, opts }.mesh(); + + EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isTriangle)); + EXPECT_EQ(24, countMeshElementsIf(staircasedMesh, isQuad)); + EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isTetrahedron)); + EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isHexahedron)); + +} + +TEST_F(StaircaseMesherTest, meshesSelectedNonzeroVolumeGroupWithHexahedra) +{ + Mesh mesh = buildCubeSurfaceMesh(0.5); + mesh.groups[0].name = "surface"; + mesh.groups.push_back(mesh.groups[0]); + mesh.groups[1].name = "volume"; + StaircaseMesherOptions options; + options.volumeGroups.insert(1); + + const Mesh result = StaircaseMesher(mesh, 4, options).mesh(); + + EXPECT_EQ("surface", result.groups[0].name); + EXPECT_EQ("volume", result.groups[1].name); + EXPECT_EQ(24, std::count_if( + result.groups[0].elements.begin(), result.groups[0].elements.end(), isQuad)); + EXPECT_EQ(4, std::count_if( + result.groups[1].elements.begin(), result.groups[1].elements.end(), isHexahedron)); +} + #if APP_LOADED -TEST_F(StaircaseMesherTest, fills_closed_volume_with_quads) +TEST_F(StaircaseMesherTest, fillsSphereAsSingleClosedUnitHexahedralVolume) { auto mesh = vtkIO::readInputMesh("testData/cases/sphere/sphere.stl"); @@ -379,20 +584,40 @@ TEST_F(StaircaseMesherTest, fills_closed_volume_with_quads) mesh.grid[Y] = utils::GridTools::linspace(-100.0, 100.0, 51); mesh.grid[Z] = utils::GridTools::linspace(-100.0, 100.0, 51); - meshlib::meshers::StaircaseMesherOptions opts; - opts.isVolume = false; - auto staircasedMesh = StaircaseMesher{mesh, 4, opts }.mesh(); - - opts.isVolume = true; - auto staircasedMeshVolume = StaircaseMesher{mesh, 4, opts }.mesh(); - - EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isTriangle)); - EXPECT_EQ(0, countMeshElementsIf(staircasedMesh, isTetrahedron)); + ASSERT_TRUE(isSingleClosedSurface(mesh.groups[0].elements)); - EXPECT_EQ(0, countMeshElementsIf(staircasedMeshVolume, isTriangle)); - EXPECT_EQ(0, countMeshElementsIf(staircasedMeshVolume, isTetrahedron)); + StaircaseMesherOptions options; + options.volumeGroups.insert(0); + options.splitHexahedra = true; - EXPECT_TRUE(countMeshElementsIf(staircasedMeshVolume, isQuad) > countMeshElementsIf(staircasedMesh, isQuad)); + // vtkIO::exportMeshToVTU("testData/cases/sphere/sphere.volume.before.vtk", mesh); + const Mesh result = StaircaseMesher{mesh, 4, options}.mesh(); + // vtkIO::exportMeshToVTU("testData/cases/sphere/sphere.volume.after.vtk", result); + + EXPECT_EQ(7967, countMeshElementsIf(result, isHexahedron)); + EXPECT_EQ(result.countElems(), countMeshElementsIf(result, isHexahedron)); + EXPECT_TRUE(isSingleClosedHexahedralVolume(result)); +} + +TEST_F(StaircaseMesherTest, fillsAlhambraAsSingleClosedUnitHexahedralVolume) +{ + auto mesh = vtkIO::readInputMesh("testData/cases/alhambra/alhambra.stl"); + mesh.grid[X] = utils::GridTools::linspace(-60.0, 60.0, 61); + mesh.grid[Y] = utils::GridTools::linspace(-60.0, 60.0, 61); + mesh.grid[Z] = utils::GridTools::linspace(-1.872734, 11.236404, 8); + ASSERT_TRUE(isSingleClosedSurface(mesh.groups[0].elements)); + + StaircaseMesherOptions options; + options.volumeGroups.insert(0); + options.splitHexahedra = true; + + // vtkIO::exportMeshToVTU("testData/cases/alhambra/alhambra.volume.before.vtk", mesh); + const Mesh result = StaircaseMesher{mesh, 4, options}.mesh(); + // vtkIO::exportMeshToVTU("testData/cases/alhambra/alhambra.volume.after.vtk", result); + + EXPECT_EQ(7255, countMeshElementsIf(result, isHexahedron)); + EXPECT_EQ(result.countElems(), countMeshElementsIf(result, isHexahedron)); + EXPECT_TRUE(isSingleClosedHexahedralVolume(result)); } @@ -582,4 +807,3 @@ TEST_F(StaircaseMesherTest, staircaser_reads_wires_correctly) #endif } - diff --git a/test/types/MeshTest.cpp b/test/types/MeshTest.cpp index a0e799d..2469288 100644 --- a/test/types/MeshTest.cpp +++ b/test/types/MeshTest.cpp @@ -1,7 +1,7 @@ #include "MeshTest.h" #ifdef TESSELLATOR_BOOST -TEST_F(DMesheRTypesMeshTest, serialization_deserialization) { +TEST_F(MeshTest, serialization_deserialization) { Mesh get = buildMesh(); const char* filename = "serialization_deserialization.txt"; @@ -23,4 +23,15 @@ TEST_F(DMesheRTypesMeshTest, serialization_deserialization) { } #endif +TEST_F(MeshTest, identifiesHexahedron) { + Element hexahedron( + {0, 1, 2, 3, 4, 5, 6, 7}, + Element::Type::Volume); + Element surface( + {0, 1, 2, 3, 4, 5, 6, 7}, + Element::Type::Surface); + + EXPECT_TRUE(hexahedron.isHexahedron()); + EXPECT_FALSE(surface.isHexahedron()); +} diff --git a/test/types/MeshTest.h b/test/types/MeshTest.h index 7887de3..3e6ab68 100644 --- a/test/types/MeshTest.h +++ b/test/types/MeshTest.h @@ -9,7 +9,7 @@ using namespace meshlib; -class DMesheRTypesMeshTest : public ::testing::Test { +class MeshTest : public ::testing::Test { protected: static Mesh buildMesh() { Grid grid;