diff --git a/.devcontainer/DOCKER_README.md b/.devcontainer/DOCKER_README.md new file mode 100644 index 0000000..7da37e8 --- /dev/null +++ b/.devcontainer/DOCKER_README.md @@ -0,0 +1,52 @@ +# Dev Container Setup + +This directory contains the configuration for a Docker-based development environment for Tessellator. + +## Quick Start + +### Prerequisites +- Docker installed and running +- VS Code with the "Dev Containers" extension + +### Opening the Project in Dev Container + +1. Open this repository in VS Code +2. Click the green icon in the bottom-left corner and select **"Reopen in Container"** +3. VS Code will build the Docker image and start the container automatically + +The entire workspace will be mounted at `/workspace` inside the container. + +## What's Included + +The dev container provides: +- Pre-configured C++ development environment +- CMake and Ninja build tools +- All project dependencies (VTK, Boost, CGAL via vcpkg) +- VS Code extensions: + - **C/C++ Tools** – Code navigation, IntelliSense, debugging + - **CMake Tools** – CMake project management + - **LLDB** – Debugger for C++ + +## Building and Testing + +Once inside the container: + +```bash +# Configure the project +cmake --preset docker -S . -B build + +# Build +cmake --build build -j + +# Run tests +build/bin/tessellator_tests + +# Run the application +build/bin/tessellator -i +``` + +## Notes + +- The container runs as root user +- CMake is pre-configured to skip auto-configure on file open for better performance +- Debugging with LLDB is enabled with appropriate capabilities diff --git a/.devcontainer/Dockerfile.base b/.devcontainer/Dockerfile.base new file mode 100644 index 0000000..7248eb0 --- /dev/null +++ b/.devcontainer/Dockerfile.base @@ -0,0 +1,34 @@ +FROM ubuntu:26.04 + +RUN apt-get update && apt-get install -y \ + build-essential \ + cmake \ + ninja-build \ + git \ + pkg-config \ + curl \ + wget \ + zip \ + unzip + +RUN apt-get install -y \ + libboost-graph-dev \ + libboost-program-options-dev \ + libvtk9-dev \ + nlohmann-json3-dev \ + libgtest-dev \ + libgmock-dev \ + libgmp-dev \ + libmpfr-dev + +RUN apt-get install -y \ + libeigen3-dev \ + libcgal-dev + +RUN apt-get install -y \ + clangd + +RUN chown -R ubuntu:ubuntu /home/ubuntu + +USER ubuntu +WORKDIR /home/ubuntu/dev diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json new file mode 100644 index 0000000..7847357 --- /dev/null +++ b/.devcontainer/devcontainer.json @@ -0,0 +1,28 @@ +{ + "name": "Tessellator", + "dockerFile": "Dockerfile.base", + "context": "..", + "remoteUser": "ubuntu", + "workspaceFolder": "/home/ubuntu/dev", + "mounts": [ + "source=${localWorkspaceFolder},target=/home/ubuntu/dev,type=bind" + ], + "customizations": { + "vscode": { + "extensions": [ + "ms-vscode.cpptools", + "ms-vscode.cmake-tools", + "vadimcn.vscode-lldb", + "matepek.vscode-catch2-test-adapter", + "llvm-vs-code-extensions.vscode-clangd" + ], + "settings": { + "cmake.configureOnOpen": false, + "cmake.configureOnEdit": false, + "cmake.autoSelectActiveFolder": false, + "cmake.generator": "Ninja" + } + } + }, + "runArgs": ["--cap-add=SYS_PTRACE", "--security-opt", "seccomp=unconfined"] +} diff --git a/.github/workflows/build-and-release.yml b/.github/workflows/build-and-release.yml index e7719cf..7c29fb8 100644 --- a/.github/workflows/build-and-release.yml +++ b/.github/workflows/build-and-release.yml @@ -5,7 +5,8 @@ on: branches: - main -env: +env: + VCPKG_BINARY_SOURCES: "clear;x-gha,readwrite" WINDOWS_FILENAME: opensemba-tessellator-windows-x64.tar.gz LINUX_FILENAME: opensemba-tessellator-linux.tar.gz @@ -27,7 +28,7 @@ jobs: strategy: matrix: preset: [ - {"os": windows-latest, "name": "msbuild", "filename": "windows-x64"}, + {"os": windows-2022, "name": "msbuild", "filename": "windows-x64"}, {"os": ubuntu-latest, "name": "gnu", "filename": "linux"} ] build-type: ["Release"] @@ -59,8 +60,17 @@ jobs: sudo apt-get update sudo apt-get install -y libvtk9-dev + - name: Export GitHub Actions cache environment variables + uses: actions/github-script@v7 + with: + script: | + core.exportVariable('ACTIONS_CACHE_URL', process.env.ACTIONS_CACHE_URL || ''); + core.exportVariable('ACTIONS_RUNTIME_TOKEN', process.env.ACTIONS_RUNTIME_TOKEN || ''); + - name: Setup vcpkg uses: lukka/run-vcpkg@v11 + with: + vcpkgGitCommitId: eed289f6e06a5e7a5c9e6a729671b0a56af7dd69 - name: Windows configure and build if: matrix.preset.name=='msbuild' diff --git a/.github/workflows/build-and-test.yml b/.github/workflows/build-and-test.yml index 73d442b..c6b8acd 100644 --- a/.github/workflows/build-and-test.yml +++ b/.github/workflows/build-and-test.yml @@ -5,25 +5,31 @@ on: branches: - main - dev - + concurrency: group: ${{ github.workflow }}-${{ github.event.pull_request.number || github.ref }} cancel-in-progress: true + +env: + VCPKG_BINARY_SOURCES: "clear;x-gha,readwrite" jobs: builds-and-tests: strategy: matrix: preset: [ - {"os": windows-latest, "name": "msbuild"}, - {"os": ubuntu-latest, "name": "gnu"} + {"os": windows-2022, "name": "msbuild"}, + {"os": ubuntu-latest, "name": "gnu"}, + {"os": ubuntu-latest, "name": "gnu-cgal"} ] - build-type: ["Debug", "Release"] + build-type: ["Release"] fail-fast: false name: ${{ matrix.preset.os }} / ${{matrix.preset.name}} / ${{matrix.build-type}} runs-on: ${{ matrix.preset.os }} + env: + VCPKG_BINARY_SOURCES: clear steps: - name: checkout repository @@ -45,13 +51,34 @@ jobs: sudo apt-get update sudo apt-get install -y libvtk9-dev + - name: Export GitHub Actions cache environment variables + uses: actions/github-script@v7 + with: + script: | + core.exportVariable('ACTIONS_CACHE_URL', process.env.ACTIONS_CACHE_URL || ''); + core.exportVariable('ACTIONS_RUNTIME_TOKEN', process.env.ACTIONS_RUNTIME_TOKEN || ''); + - name: Setup vcpkg uses: lukka/run-vcpkg@v11 + with: + vcpkgGitCommitId: eed289f6e06a5e7a5c9e6a729671b0a56af7dd69 - name: Windows configure and build if: matrix.preset.name=='msbuild' run: | - cmake --preset ${{matrix.preset.name}} -S . -B build + $configured = $false + for ($attempt = 1; $attempt -le 3; $attempt++) { + cmake --preset ${{matrix.preset.name}} -S . -B build + if ($LASTEXITCODE -eq 0) { + $configured = $true + break + } + if ($attempt -lt 3) { + Write-Host "Configure failed (attempt $attempt/3). Retrying in 20 seconds..." + Start-Sleep -Seconds 20 + } + } + if (-not $configured) { exit 1 } cmake --build build --config ${{matrix.build-type}} -j - name: Windows Run tests @@ -59,14 +86,19 @@ 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: | - cmake --preset ${{matrix.preset.name}} -S . -B build + configured=0 + for attempt in 1 2 3; do + cmake --preset ${{matrix.preset.name}} -S . -B build && configured=1 && break + if [ "$attempt" -lt 3 ]; then + echo "Configure failed (attempt $attempt/3). Retrying in 20 seconds..." + sleep 20 + fi + done + [ "$configured" -eq 1 ] || exit 1 cmake --build build -j - name: Ubuntu Run tests - if: matrix.preset.name=='gnu' + if: matrix.preset.os=='ubuntu-latest' run: build/bin/tessellator_tests - - - \ No newline at end of file diff --git a/.gitignore b/.gitignore index 2c7a665..130eeec 100644 --- a/.gitignore +++ b/.gitignore @@ -8,11 +8,15 @@ build/ .settings .project .cproject +.cache src/*.json .vs/ -.vscode/ +.vscode/settings.json +.vscode/settings.dev.json +.vscode/launch.json +.vscode/tasks.json testData/*_out.stl *_out.stl @@ -25,4 +29,5 @@ CMakeUserPresets.json sliced.vtk contour.vtk -testData/ \ No newline at end of file +testData/ +build diff --git a/.vscode/launch.json b/.vscode/launch.json new file mode 100644 index 0000000..96b5b7a --- /dev/null +++ b/.vscode/launch.json @@ -0,0 +1,105 @@ +{ + "version": "0.2.0", + "configurations": [ + { + "name": "Debug Tessellator (Docker)", + "type": "lldb", + "request": "launch", + "program": "${workspaceFolder}/build-dbg/bin/tessellator", + "args": ["-i", "${workspaceFolder}/testData/cases/cone/cone.tessellator.json"], + "stopOnEntry": false, + "cwd": "${workspaceFolder}", + "env": {} + }, + { + "name": "Debug Tests (Docker)", + "type": "lldb", + "request": "launch", + "program": "${workspaceFolder}/build-dbg/bin/tessellator_tests", + "args": ["--gtest_filter=*"], + "stopOnEntry": false, + "cwd": "${workspaceFolder}", + "env": {}, + "sourceMap": {} + }, + { + "name": "Debug Tests (Docker) - Specific Suite", + "type": "lldb", + "request": "launch", + "program": "${workspaceFolder}/build-dbg/bin/tessellator_tests", + "args": ["--gtest_filter=${input:testSuite}"], + "stopOnEntry": false, + "cwd": "${workspaceFolder}", + "env": {}, + "sourceMap": {}, + }, + { + "name": "tessellator (gdb)", + "type": "cppdbg", + "request": "launch", + "program": "${workspaceFolder}/build-dbg/bin/tessellator", + "args": ["-i", "testData/cases/alhambra/alhambra.tessellator.json"], + "stopAtEntry": false, + "cwd": "${workspaceFolder}", + "environment": [], + "externalConsole": false, + "MIMode": "gdb", + "visualizerFile": [ + "${workspaceFolder}/resources/Eigen.natvis", + "${workspaceFolder}/resources/nlohmann_json.natvis" + ], + "additionalSOLibSearchPath": "", + "showDisplayString": true, + "setupCommands": [ + { + "description": "Enable pretty-printing for gdb", + "text": "-enable-pretty-printing", + "ignoreFailures": true + }, + { + "description": "Set Disassembly Flavor to Intel", + "text": "-gdb-set disassembly-flavor intel", + "ignoreFailures": true + } + ] + }, + { + "name": "tessellator_tests (gdb)", + "type": "cppdbg", + "request": "launch", + "program": "${workspaceFolder}/build-dbg/bin/tessellator_tests", + "args": [], + "stopAtEntry": false, + "cwd": "${workspaceFolder}", + "environment": [], + "externalConsole": false, + "MIMode": "gdb", + "visualizerFile": [ + "${workspaceFolder}/resources/Eigen.natvis", + "${workspaceFolder}/resources/nlohmann_json.natvis" + ], + "additionalSOLibSearchPath": "", + "showDisplayString": true, + "setupCommands": [ + { + "description": "Enable pretty-printing for gdb", + "text": "-enable-pretty-printing", + "ignoreFailures": true + }, + { + "description": "Set Disassembly Flavor to Intel", + "text": "-gdb-set disassembly-flavor intel", + "ignoreFailures": true + } + ] + } + ], + "inputs": [ + { + "id": "testSuite", + "type": "promptString", + "description": "Enter test filter (e.g., StaircaseMesherTest.* or StaircaseMesherTest.TestName)", + "default": "*" + } + ] +} \ No newline at end of file diff --git a/.vscode/settings.dev.json b/.vscode/settings.dev.json new file mode 100644 index 0000000..9d190ee --- /dev/null +++ b/.vscode/settings.dev.json @@ -0,0 +1,15 @@ +{ + "cmake.configureOnOpen": false, + "cmake.configureOnEdit": false, + "cmake.autoSelectActiveFolder": false, + "testMate.cpp.test.executables": "{build,build-dbg,Build,BUILD,out,Out,OUT}/**/*{test,Test,TEST}*", + "testMate.cpp.test.workingDirectory": "${workspaceFolder}", + "extensions.ignoreRecommendations": true, + "C_Cpp.intelliSenseEngine": "disabled", + "clangd.arguments": ["-log=verbose", + "-pretty", + "--background-index", + "--query-driver=/usr/bin/g++" + ], + "clangd.fallbackFlags": ["-std=c++17"] +} diff --git a/CITATION.cff b/CITATION.cff index 5368c9d..86152ca 100644 --- a/CITATION.cff +++ b/CITATION.cff @@ -12,7 +12,7 @@ authors: orcid: https://orcid.org/0000-0001-7317-1423 - family-names: Rubio Bretones given-names: Amelia - orcid: https://orcid.org/00000-0002-9337-9093 + orcid: https://orcid.org/0000-0002-9337-9093 url: "https://github.com/opensemba/tessellator" title: "opensemba/tessellator an FDTD meshing tool" date-released: 2024-07-10 \ No newline at end of file diff --git a/CLAUDE.md b/CLAUDE.md new file mode 100644 index 0000000..f9e2bac --- /dev/null +++ b/CLAUDE.md @@ -0,0 +1,205 @@ +# CLAUDE.md + +This file provides guidance to Claude Code (claude.ai/code) when working with code in this repository. + +## Overview + +**Tessellator** is a C++17 mesher designed to generate meshes and data structures optimized for FDTD algorithms. The project supports multiple meshing strategies (staircased, conformal, offgrid) and includes extensive test coverage using Google Test. + +## Key Architecture + +### Core Components + +The codebase is organized into modular namespaces under `src/`: + +- **`meshers/`** – Abstract mesh generation interface with three implementations: + - `StaircaseMesher` – Generates staircased meshes from geometric inputs + - `ConformalMesher` – Creates conformal meshes with fixed-distance grid plane intersections + - `OffgridMesher` – Produces offgrid meshes (graded rectilinear support) + - All inherit from `MesherBase` which handles input mesh processing and grid construction + +- **`core/`** – Fundamental mesh processing algorithms: + - `Slicer` – Slices input geometry along grid planes + - `Staircaser` – Converts sliced geometry into staircase patterns + - `Smoother` – Post-processing to smooth mesh irregularities + - `Collapser` – Merges adjacent elements for optimization + - `Snapper` – Aligns geometry with grid coordinates + +- **`types/`** – Core data structures: + - `Mesh` – Central data structure holding coordinates, elements (nodes/lines/surfaces/volumes), and metadata + - `Mesher` – Abstract interface that all meshers implement + - `Vector` – Template-based 3D vector type + - `Grid` – Rectilinear grid representation (3D array of coordinate positions) + +- **`cgal/`** – Optional CGAL-based advanced geometry operations (controlled by `TESSELLATOR_ENABLE_CGAL`): + - `Filler` – Fills closed polyhedrons using CGAL + - `Delaunator` – Triangulation operations + - `Manifolder` – Manifold repair utilities + - `Repairer` – Polyhedron repair tools + - `HPolygonSet`, `PolyhedronTools` – Polygon/polyhedron manipulation + +- **`utils/`** – Utility functions: + - `MeshTools`, `GridTools` – Mesh/grid manipulation + - `Geometry` – Geometric calculations + - `ConvexHull` – Convex hull computation + - `CoordGraph`, `ElemGraph` – Topological graph representations + - `RedundancyCleaner` – Removes redundant geometric data + +- **`app/`** – Application interface: + - `launcher` – CLI entry point and argument parsing + - `vtkIO` – VTK file import/export (STL, VTK formats) + +### Data Flow + +Input geometry (STL/VTK) → VTK I/O → Mesher (chooses algorithm) → Core processing (slicer/staircaser/smoother) → Grid alignment → Output mesh (VTK) + +## Build System + +Uses **CMake 3.20+** with presets and vcpkg for dependency management. + +### Build Presets + +```bash +# GNU/Linux (Ninja) +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 +``` + +### CMake Options + +- `TESSELLATOR_ENABLE_TESTS` (ON by default) – Build test suite +- `TESSELLATOR_ENABLE_CGAL` (OFF by default) – Enable CGAL-based geometry operations +- `TESSELLATOR_EXECUTION_POLICIES` (OFF by default) – Parallel execution policies + +### Dependencies + +Managed via vcpkg manifest: +- **Required**: VTK, Boost +- **Optional**: CGAL (if `TESSELLATOR_ENABLE_CGAL=ON`) + +To set up locally, create a `CMakeUserPreset.json` file: + +```json +{ + "version": 4, + "include": ["CMakePresets.json"], + "configurePresets": [ + { + "name": "gnu-local", + "environment": { + "VCPKG_ROOT": "~/workspace/vcpkg/" + }, + "cacheVariables": { + "TESSELLATOR_ENABLE_CGAL": false + }, + "inherits": "gnu" + } + ] +} +``` + +## Common Development Tasks + +### Build + +```bash +cmake --preset gnu -S . -B build +cmake --build build -j +``` + +Output binaries go to `build/bin/` and libraries to `build/lib/`. + +### Run Tests + +```bash +# Run all tests +build/bin/tessellator_tests + +# Run specific test suite (e.g., all StaircaseMesher tests) +build/bin/tessellator_tests --gtest_filter=StaircaseMesher* + +# Run single test method +build/bin/tessellator_tests --gtest_filter=StaircaseMesherTest.TestName +``` + +Tests are organized by module under `test/` mirroring the source structure. + +### Run Application + +```bash +# Build and run with a tessellator JSON input file +build/bin/tessellator -i mesh_definition.tessellator.json +``` + +## Code Style & Patterns + +- **C++17** standard; no C++20+ features +- **Namespaces**: Organized as `meshlib`, `meshlib::meshers`, `meshlib::core`, etc. +- **Mesh representation**: All geometry flows through the `Mesh` object with its `Element` struct (type: Node/Line/Surface/Volume) +- **Grid**: Stored as `Grid = std::array, 3>` (coordinate arrays per axis) +- **Logging**: Use static `MesherBase::log()` methods with indentation level for hierarchical output +- **Vertex references**: Elements reference coordinates by index (`CoordinateId`) + +## Testing + +Google Test (GTest) framework with fixtures. Tests for each module live in `test//`. Key fixtures: + +- `MeshFixtures.h` – Shared mesh construction utilities for tests + +**Important**: When modifying core algorithms (Slicer, Staircaser, Smoother, Snapper, Collapser), update corresponding tests in `test/core/` and `test/meshers/`. + +## JSON Input Format + +The tessellator expects a JSON file with two main entries: + +### Grid Definition + +Define the grid structure using either uniform or rectilinear (graded) cells: + +```json +"grid": { + "numberOfCells": [20, 20, 30], + "boundingBox": [[-1, -1, -1], [1, 1, 2]] +} +``` + +Or for non-uniform grids: + +```json +"grid": { + "planes": [ + [600, 603.25], + [25.0, 30.5, 92, 130, 1000], + [1000, 1111, 1111.1] + ] +} +``` + +### Object Definition + +```json +"object": {"filename": "geometry.stl"} +``` + +The filename is relative to the JSON file location. + +## File Organization + +- **Headers** (`.h`) – Include guards, template implementations +- **Implementations** (`.cpp`) – Separate compilation units +- **Test files** – Follow `Test.cpp` pattern with GTest syntax + +## Known Constraints + +- CGAL-dependent code is conditionally compiled; check `if (TESSELLATOR_ENABLE_CGAL)` guards when making changes +- Grid must be rectilinear; arbitrary unstructured grids are not supported +- Mesh must be manifold for some operations (see `Manifolder`) +- Elements reference coordinates by index, not direct storage diff --git a/CMakeLists.txt b/CMakeLists.txt index c9c568f..5d3a1b0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,26 +1,67 @@ cmake_minimum_required(VERSION 3.20) +# Toolchain detection: support both vcpkg (local dev) and system libraries (Docker) if(NOT DEFINED CMAKE_TOOLCHAIN_FILE AND DEFINED ENV{CMAKE_TOOLCHAIN_FILE}) set(CMAKE_TOOLCHAIN_FILE $ENV{CMAKE_TOOLCHAIN_FILE}) endif() +# If CMAKE_TOOLCHAIN_FILE is set, verify it exists (avoid silent failures) +if(DEFINED CMAKE_TOOLCHAIN_FILE AND NOT EXISTS "${CMAKE_TOOLCHAIN_FILE}") + message(WARNING "CMAKE_TOOLCHAIN_FILE specified but not found: ${CMAKE_TOOLCHAIN_FILE}") + message(STATUS "Falling back to system libraries (Docker mode)") + unset(CMAKE_TOOLCHAIN_FILE) +endif() + +if(DEFINED CMAKE_TOOLCHAIN_FILE) + message(STATUS "Using vcpkg toolchain: ${CMAKE_TOOLCHAIN_FILE}") +else() + message(STATUS "Using system libraries (vcpkg toolchain not detected)") +endif() + set(CMAKE_CXX_STANDARD 17) set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib) set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib) set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin) +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) +if(DOCKER_EXPORT_COMPILE_COMMANDS) + set(CMAKE_EXPORT_COMPILE_COMMANDS ON) +endif() if(TESSELLATOR_ENABLE_CGAL) list(APPEND VCPKG_MANIFEST_FEATURES "cgal") + SET( EIGEN3_INCLUDE_DIR "$ENV{EIGEN3_INCLUDE_DIR}" ) + if (EIGEN3_INCLUDE_DIR) + INCLUDE_DIRECTORIES ( "$ENV{EIGEN3_INCLUDE_DIR}" ) + INCLUDE_DIRECTORIES ( "$ENV{EIGEN3_INCLUDE_DIR}/Eigen" ) + endif() endif() project(tessellator CXX) +find_package(VTK COMPONENTS + CommonCore + IOGeometry + IOLegacy + FiltersCore +) + +if(VTK_FOUND) + set(TESSELLATOR_LOAD_APP ON) + add_definitions(-DAPP_LOADED=1) + +else() + message(STATUS "VTK not found - tessellator app will not be built") + add_definitions(-DAPP_LOADED=0) +endif() + add_subdirectory(src/) if(TESSELLATOR_ENABLE_TESTS) diff --git a/CMakePresets.json b/CMakePresets.json index d80f52a..93003ba 100644 --- a/CMakePresets.json +++ b/CMakePresets.json @@ -1,38 +1,106 @@ { - "version": 4, - "configurePresets": [ - { - "name": "default", - "hidden": true, - "binaryDir": "build/", - "cacheVariables": { - "CMAKE_TOOLCHAIN_FILE": { - "type": "FILEPATH", - "value": "$env{VCPKG_ROOT}/scripts/buildsystems/vcpkg.cmake" + "version": 4, + "configurePresets": [ + { + "name": "default", + "hidden": true, + "binaryDir": "build/", + "cacheVariables": { + "CMAKE_TOOLCHAIN_FILE": { + "type": "FILEPATH", + "value": "$env{VCPKG_ROOT}/scripts/buildsystems/vcpkg.cmake" + }, + "TESSELLATOR_ENABLE_CGAL": false + } + }, + { + "name": "msbuild", + "displayName": "MSBuild Configure Settings", + "generator": "Visual Studio 17 2022", + "architecture": { + "strategy": "set", + "value": "x64" + }, + "inherits": "default" + }, + { + "name": "gnu", + "displayName": "GNU g++ compiler", + "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", + "inherits": "gnu", + "binaryDir": "build-dbg/", + "cacheVariables": { + "CMAKE_BUILD_TYPE": "Debug" + } + }, + { + "name": "docker", + "displayName": "Docker (system libraries)", + "description": "Uses system-installed libraries, skips vcpkg. Optimized for containerized builds.", + "generator": "Ninja", + "binaryDir": "build/", + "cacheVariables": { + "CMAKE_CXX_COMPILER": "g++", + "CMAKE_PREFIX_PATH": "/usr/local", + "CMAKE_FIND_ROOT_PATH": "/usr/local", + "TESSELLATOR_ENABLE_TESTS": "ON", + "TESSELLATOR_ENABLE_CGAL": "OFF", + "DOCKER_EXPORT_COMPILE_COMMANDS": "ON" + }, + "environment": { + "EIGEN3_INCLUDE_DIR": "/usr/include/eigen3" + } + }, + { + "name": "docker-dbg", + "inherits": "docker", + "displayName": "Docker-dbg (system libraries)", + "description": "Uses system-installed libraries, skips vcpkg. Optimized for containerized builds. Debug mode", + "binaryDir": "build-dbg/", + "cacheVariables": { + "CMAKE_BUILD_TYPE": "Debug", + "CMAKE_CXX_COMPILER": "g++", + "CMAKE_PREFIX_PATH": "/usr/local", + "CMAKE_FIND_ROOT_PATH": "/usr/local", + "TESSELLATOR_ENABLE_TESTS": "ON", + "TESSELLATOR_ENABLE_CGAL": "OFF", + "DOCKER_EXPORT_COMPILE_COMMANDS": "ON", + "CMAKE_CXX_FLAGS_INIT": "-g3" + }, + "environment": { + "EIGEN3_INCLUDE_DIR": "/usr/include/eigen3" + } } - } - }, - { - "name": "msbuild", - "displayName": "MSBuild Configure Settings", - "generator": "Visual Studio 17 2022", - "architecture": { - "strategy": "set", - "value": "x64" - }, - "inherits": "default" - }, - { - "name": "gnu", - "displayName": "GNU g++ compiler", - "generator": "Ninja", - "inherits": "default" - } - ], - "buildPresets": [ - { - "name": "default", - "configurePreset": "default" - } - ] -} \ No newline at end of file + ], + "buildPresets": [ + { + "name": "default", + "configurePreset": "default" + }, + { + "name": "docker", + "configurePreset": "docker" + }, + { + "name": "docker-dbg", + "configurePreset": "docker-dbg" + }, + { + "name": "gnu-cgal", + "configurePreset": "gnu-cgal" + } + ] +} diff --git a/README.md b/README.md index b777cec..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" } @@ -50,13 +51,13 @@ The main binary is `tessellator`, which uses a tessellator json format, which wi ``` ## JSON Format -The two main entries are as follows: +The main entries are as follows: ### `` This object must always be present and contains the structure of the grid, which will be used to slice and adjust the mesh provided. It must contain one of these two sets of entries: - ``: is an array of three positive integers which indicate the number of cells in each Cartesian direction. In case of having this entry, it also must contain a ``: - - `` is represented by an array which contairs two triplets of integers, representing the minimum and maximum values of the gread in each cartesian direction. + - `` is represented by an array which contains two triplets of integers, representing the minimum and maximum values of the grid in each cartesian direction. ```json "grid": { @@ -80,17 +81,101 @@ This object must always be present and contains the structure of the grid, which } ``` -### `` -This contains the information about the mesh file. It must contain the following entry: +### `` or `` +This contains the information about the mesh file(s). You can specify a single object or multiple objects: -- `filename`: with an string containing the name of the mesh file. Its location is relative to that of the json file. - - Example: +**Single object:** +- `filename`: A string containing the name of the mesh file. Its location is relative to that of the json file. ```json "object": {"filename": "thinCylinder.stl"} ``` +**Multiple objects:** +- `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) + - `mesher`: (optional) Override the global mesher settings for this specific object + +```json + "objects": [ + {"filename": "object1.stl", "group": "group1"}, + {"filename": "object2.stl", "group": "group2", "mesher": {"type": "conformal"}} + ] +``` + +### `` +This optional entry configures the meshing algorithm and its options. If not specified, the staircase mesher is used with default options. + +**Mesher types:** +- `staircase` (default): Generates staircased meshes from geometric inputs +- `conformal`: Creates conformal meshes with fixed-distance grid plane intersections + +**Mesher options:** + +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 +- `forbiddenLength`: Minimum length threshold for snapping + +**Global options:** +- `exportGrid`: (boolean, default: true) Controls whether to export the grid file + +Example with staircase mesher and compression enabled: +```json + "mesher": { + "type": "staircase", + "options": { + "compress": true, + "exportGrid": true + } + } +``` + +Example with conformal mesher: +```json + "mesher": { + "type": "conformal", + "options": { + "edgePoints": true, + "forbiddenLength": 0.001 + } + } +``` + +### Output Files +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.grid.vtk` - Grid file (if `exportGrid` is true) + +### Complete Example +```json +{ + "grid": { + "numberOfCells": [50, 50, 50], + "boundingBox": [ + [-100.0, -100.0, -100.0], + [ 100.0, 100.0, 100.0] + ] + }, + "objects": [ + {"filename": "sphere.stl"}, + {"filename": "cylinder.stl", "mesher": {"type": "conformal"}} + ], + "mesher": { + "type": "staircase", + "options": { + "compress": true, + "exportGrid": true + } + } +} +``` + ## Contributing ## Citing this work diff --git a/overlayPorts/libaec/README.md b/overlayPorts/libaec/README.md new file mode 100644 index 0000000..ca36a07 --- /dev/null +++ b/overlayPorts/libaec/README.md @@ -0,0 +1,12 @@ +Original port: https://github.com/microsoft/vcpkg/tree/fedfb6ae868887a13f8a0ba3b29603bd7eb9f118/ports/libaec + +A different download url is used to avoid connection errors. + +libaec provides CMake targets: +``` Cmake + find_package(libaec CONFIG REQUIRED) + # libaec API + target_link_libraries(main PRIVATE libaec::aec) + # szip compatible API + target_link_libraries(main PRIVATE libaec::sz) +``` \ No newline at end of file diff --git a/overlayPorts/libaec/portfile.cmake b/overlayPorts/libaec/portfile.cmake new file mode 100644 index 0000000..cdd5c9f --- /dev/null +++ b/overlayPorts/libaec/portfile.cmake @@ -0,0 +1,44 @@ +vcpkg_from_github( + OUT_SOURCE_PATH SOURCE_PATH + REPO Deutsches-Klimarechenzentrum/libaec + REF v${VERSION} + SHA512 76df7501d1b7d91a43b525ba828f092f18d83f8ab09a9331e5758f93942a9758ad580baca8f9316b92a98639bde2e23cacbc2f33f52d0dd98ce7efe412cf43cd + HEAD_REF master +) + +string(COMPARE EQUAL "${VCPKG_LIBRARY_LINKAGE}" "static" BUILD_STATIC) + +vcpkg_cmake_configure( + SOURCE_PATH "${SOURCE_PATH}" + OPTIONS + -DBUILD_STATIC_LIBS=${BUILD_STATIC} + -Dlibaec_INSTALL_CMAKEDIR=share/${PORT} +) +vcpkg_cmake_install() +vcpkg_copy_pdbs() +vcpkg_cmake_config_fixup() +vcpkg_replace_string("${CURRENT_PACKAGES_DIR}/share/libaec/libaec-config.cmake" + "if(libaec_USE_STATIC_LIBS)" + "if(\"${BUILD_STATIC}\") # forced by vcpkg" +) + +# Compatibility with user's CMake < 3.18 (vcpkg claims support for >= 3.16): +# Make imported targets global so that libaec-config.cmake can create ALIAS targets. +set(_target_file "libaec_shared-targets") +if(BUILD_STATIC) + set(_target_file "libaec_static-targets") +endif() +file(READ "${CURRENT_PACKAGES_DIR}/share/libaec/${_target_file}.cmake" libaec_targets) +string(REGEX REPLACE " (SHARED|STATIC) IMPORTED" " \\1 IMPORTED \${libaec_maybe_global}" libaec_targets "${libaec_targets}") +file(WRITE "${CURRENT_PACKAGES_DIR}/share/libaec/${_target_file}.cmake" "set(libaec_maybe_global \"\") +if(CMAKE_VERSION VERSION_LESS 3.18) + set(libaec_maybe_global \"GLOBAL\") +endif() +${libaec_targets} +" +) + +file(REMOVE_RECURSE "${CURRENT_PACKAGES_DIR}/debug/include") + +file(INSTALL "${CURRENT_PORT_DIR}/usage" DESTINATION "${CURRENT_PACKAGES_DIR}/share/${PORT}") +vcpkg_install_copyright(FILE_LIST "${SOURCE_PATH}/LICENSE.txt") \ No newline at end of file diff --git a/overlayPorts/libaec/usage b/overlayPorts/libaec/usage new file mode 100644 index 0000000..26ebd09 --- /dev/null +++ b/overlayPorts/libaec/usage @@ -0,0 +1,7 @@ +libaec provides CMake targets: + + find_package(libaec CONFIG REQUIRED) + # libaec API + target_link_libraries(main PRIVATE libaec::aec) + # szip compatible API + target_link_libraries(main PRIVATE libaec::sz) \ No newline at end of file diff --git a/overlayPorts/libaec/vcpkg.json b/overlayPorts/libaec/vcpkg.json new file mode 100644 index 0000000..a4f3d55 --- /dev/null +++ b/overlayPorts/libaec/vcpkg.json @@ -0,0 +1,17 @@ +{ + "name": "libaec", + "version": "1.1.6", + "description": "Adaptive Entropy Coding library", + "homepage": "https://gitlab.dkrz.de/k202009/libaec", + "license": "BSD-2-Clause", + "dependencies": [ + { + "name": "vcpkg-cmake", + "host": true + }, + { + "name": "vcpkg-cmake-config", + "host": true + } + ] +} \ No newline at end of file diff --git a/resources/Eigen.natvis b/resources/Eigen.natvis new file mode 100644 index 0000000..dbed89e --- /dev/null +++ b/resources/Eigen.natvis @@ -0,0 +1,253 @@ + + + + + + + + [{$T2}, {$T3}] + + + 2 + $i==0 ? $T2 : $T3 + m_storage.m_data.array + + + Backward + 2 + $i==0 ? $T2 : $T3 + m_storage.m_data.array + + + + + + + [2, 2] + + + {m_storage.m_data.array[0]} {m_storage.m_data.array[1]} + + + {m_storage.m_data.array[0]} {m_storage.m_data.array[2]} + + + {m_storage.m_data.array[2]} {m_storage.m_data.array[3]} + + + {m_storage.m_data.array[1]} {m_storage.m_data.array[3]} + + + + + + + [3, 3] + + + {m_storage.m_data.array[0]} {m_storage.m_data.array[1]} {m_storage.m_data.array[2]} + + + {m_storage.m_data.array[0]} {m_storage.m_data.array[3]} {m_storage.m_data.array[6]} + + + {m_storage.m_data.array[3]} {m_storage.m_data.array[4]} {m_storage.m_data.array[5]} + + + {m_storage.m_data.array[1]} {m_storage.m_data.array[4]} {m_storage.m_data.array[7]} + + + {m_storage.m_data.array[6]} {m_storage.m_data.array[7]} {m_storage.m_data.array[8]} + + + {m_storage.m_data.array[2]} {m_storage.m_data.array[5]} {m_storage.m_data.array[8]} + + + + + + + [3, 4] + + + {m_storage.m_data.array[0]} {m_storage.m_data.array[1]} {m_storage.m_data.array[2]} {m_storage.m_data.array[3]} + + + {m_storage.m_data.array[0]} {m_storage.m_data.array[3]} {m_storage.m_data.array[6]} {m_storage.m_data.array[9]} + + + {m_storage.m_data.array[4]} {m_storage.m_data.array[5]} {m_storage.m_data.array[6]} {m_storage.m_data.array[7]} + + + {m_storage.m_data.array[1]} {m_storage.m_data.array[4]} {m_storage.m_data.array[7]} {m_storage.m_data.array[10]} + + + {m_storage.m_data.array[8]} {m_storage.m_data.array[9]} {m_storage.m_data.array[10]} {m_storage.m_data.array[11]} + + + {m_storage.m_data.array[2]} {m_storage.m_data.array[5]} {m_storage.m_data.array[8]} {m_storage.m_data.array[11]} + + + + + + + [4, 4] + + + {m_storage.m_data.array[0]} {m_storage.m_data.array[1]} {m_storage.m_data.array[2]} {m_storage.m_data.array[3]} + + + {m_storage.m_data.array[0]} {m_storage.m_data.array[4]} {m_storage.m_data.array[8]} {m_storage.m_data.array[12]} + + + {m_storage.m_data.array[4]} {m_storage.m_data.array[5]} {m_storage.m_data.array[6]} {m_storage.m_data.array[7]} + + + {m_storage.m_data.array[1]} {m_storage.m_data.array[5]} {m_storage.m_data.array[9]} {m_storage.m_data.array[13]} + + + {m_storage.m_data.array[8]} {m_storage.m_data.array[9]} {m_storage.m_data.array[10]} {m_storage.m_data.array[11]} + + + {m_storage.m_data.array[2]} {m_storage.m_data.array[6]} {m_storage.m_data.array[10]} {m_storage.m_data.array[14]} + + + {m_storage.m_data.array[12]} {m_storage.m_data.array[13]} {m_storage.m_data.array[14]} {m_storage.m_data.array[15]} + + + {m_storage.m_data.array[3]} {m_storage.m_data.array[7]} {m_storage.m_data.array[11]} {m_storage.m_data.array[15]} + + + + + + + + empty + [{m_storage.m_rows}, {m_storage.m_cols}] (dynamic matrix) + + + 2 + $i==0 ? m_storage.m_rows : m_storage.m_cols + m_storage.m_data + + + Backward + 2 + $i==0 ? m_storage.m_rows : m_storage.m_cols + m_storage.m_data + + + + + + + empty + [{$T2}, {m_storage.m_cols}] (dynamic column matrix) + + + 2 + $i==0 ? $T2 : m_storage.m_cols + m_storage.m_data + + + Backward + 2 + $i==0 ? $T2 : m_storage.m_cols + m_storage.m_data + + + + + + + + empty + [{m_storage.m_rows}, {$T2}] (dynamic row matrix) + + + 2 + $i==0 ? m_storage.m_rows : $T2 + m_storage.m_data + + + Backward + 2 + $i==0 ? m_storage.m_rows : $T2 + m_storage.m_data + + + + + + + + empty + [{m_storage.m_cols}] (dynamic column vector) + + m_storage.m_cols + + m_storage.m_cols + m_storage.m_data + + + + + + + + empty + [{m_storage.m_rows}] (dynamic row vector) + + m_storage.m_rows + + m_storage.m_rows + m_storage.m_data + + + + + + + + [1] {m_storage.m_data.array[0]} + + m_storage.m_data.array[0] + + + + + + + [2] {m_storage.m_data.array[0]} {m_storage.m_data.array[1]} + + m_storage.m_data.array[0] + m_storage.m_data.array[1] + + + + + + + [3] {m_storage.m_data.array[0]} {m_storage.m_data.array[1]} {m_storage.m_data.array[2]} + + m_storage.m_data.array[0] + m_storage.m_data.array[1] + m_storage.m_data.array[2] + + + + + + + [4] {m_storage.m_data.array[0]} {m_storage.m_data.array[1]} {m_storage.m_data.array[2]} {m_storage.m_data.array[3]} + + m_storage.m_data.array[0] + m_storage.m_data.array[1] + m_storage.m_data.array[2] + m_storage.m_data.array[3] + + + + diff --git a/resources/nlohmann_json.natvis b/resources/nlohmann_json.natvis new file mode 100644 index 0000000..5449cae --- /dev/null +++ b/resources/nlohmann_json.natvis @@ -0,0 +1,278 @@ + + + + + + + + + + null + {*(m_value.object)} + {*(m_value.array)} + {*(m_value.string)} + {m_value.boolean} + {m_value.number_integer} + {m_value.number_unsigned} + {m_value.number_float} + discarded + + + *(m_value.object),view(simple) + + + *(m_value.array),view(simple) + + + + + + + {second} + + second + + + + + + null + {*(m_value.object)} + {*(m_value.array)} + {*(m_value.string)} + {m_value.boolean} + {m_value.number_integer} + {m_value.number_unsigned} + {m_value.number_float} + discarded + + + *(m_value.object),view(simple) + + + *(m_value.array),view(simple) + + + + + + + {second} + + second + + + + + + null + {*(m_value.object)} + {*(m_value.array)} + {*(m_value.string)} + {m_value.boolean} + {m_value.number_integer} + {m_value.number_unsigned} + {m_value.number_float} + discarded + + + *(m_value.object),view(simple) + + + *(m_value.array),view(simple) + + + + + + + {second} + + second + + + + + + null + {*(m_value.object)} + {*(m_value.array)} + {*(m_value.string)} + {m_value.boolean} + {m_value.number_integer} + {m_value.number_unsigned} + {m_value.number_float} + discarded + + + *(m_value.object),view(simple) + + + *(m_value.array),view(simple) + + + + + + + {second} + + second + + + + + + null + {*(m_value.object)} + {*(m_value.array)} + {*(m_value.string)} + {m_value.boolean} + {m_value.number_integer} + {m_value.number_unsigned} + {m_value.number_float} + discarded + + + *(m_value.object),view(simple) + + + *(m_value.array),view(simple) + + + + + + + {second} + + second + + + + + + null + {*(m_value.object)} + {*(m_value.array)} + {*(m_value.string)} + {m_value.boolean} + {m_value.number_integer} + {m_value.number_unsigned} + {m_value.number_float} + discarded + + + *(m_value.object),view(simple) + + + *(m_value.array),view(simple) + + + + + + + {second} + + second + + + + + + null + {*(m_value.object)} + {*(m_value.array)} + {*(m_value.string)} + {m_value.boolean} + {m_value.number_integer} + {m_value.number_unsigned} + {m_value.number_float} + discarded + + + *(m_value.object),view(simple) + + + *(m_value.array),view(simple) + + + + + + + {second} + + second + + + + + + null + {*(m_value.object)} + {*(m_value.array)} + {*(m_value.string)} + {m_value.boolean} + {m_value.number_integer} + {m_value.number_unsigned} + {m_value.number_float} + discarded + + + *(m_value.object),view(simple) + + + *(m_value.array),view(simple) + + + + + + + {second} + + second + + + + + + null + {*(m_value.object)} + {*(m_value.array)} + {*(m_value.string)} + {m_value.boolean} + {m_value.number_integer} + {m_value.number_unsigned} + {m_value.number_float} + discarded + + + *(m_value.object),view(simple) + + + *(m_value.array),view(simple) + + + + + + + {second} + + second + + + + diff --git a/src/app/CMakeLists.txt b/src/app/CMakeLists.txt index 6012a8c..ed3d2c7 100644 --- a/src/app/CMakeLists.txt +++ b/src/app/CMakeLists.txt @@ -1,29 +1,33 @@ message(STATUS "Creating build system for tessellator-app") -add_library(tessellator-app - "vtkIO.cpp" - "launcher.cpp" -) +if(TESSELLATOR_LOAD_APP) + find_package(VTK COMPONENTS + CommonCore + IOGeometry + IOLegacy + FiltersCore + ) -find_package(VTK COMPONENTS - CommonCore - IOGeometry - IOLegacy - FiltersCore -) + add_library(tessellator-app + "vtkIO.cpp" + "launcher.cpp" + ) -find_package(Boost COMPONENTS program_options) + find_package(Boost COMPONENTS program_options) -find_package(nlohmann_json) + find_package(nlohmann_json) -target_link_libraries(tessellator-app - ${VTK_LIBRARIES} - Boost::program_options - nlohmann_json::nlohmann_json -) + target_link_libraries(tessellator-app + ${VTK_LIBRARIES} + Boost::program_options + nlohmann_json::nlohmann_json + ) + add_executable(tessellator + "tessellator.cpp" + ) -add_executable(tessellator - "tessellator.cpp" -) - -target_link_libraries(tessellator tessellator-app tessellator-meshers) + target_link_libraries(tessellator tessellator-app tessellator-meshers) + + else() + add_library(tessellator-app INTERFACE) +endif() diff --git a/src/app/launcher.cpp b/src/app/launcher.cpp index b8dd444..5221e49 100644 --- a/src/app/launcher.cpp +++ b/src/app/launcher.cpp @@ -1,9 +1,11 @@ #include "launcher.h" #include "vtkIO.h" +#include "meshers/MesherBase.h" #include "meshers/StaircaseMesher.h" #include "meshers/ConformalMesher.h" #include "utils/GridTools.h" +#include "utils/MeshTools.h" #include #include @@ -13,6 +15,7 @@ #include #include #include +#include namespace meshlib::app { @@ -21,20 +24,20 @@ using namespace vtkIO; namespace po = boost::program_options; -Grid parseGridFromJSON(const nlohmann::json &j) +Grid parseGridFromJSON(const nlohmann::json &fileData) { - if (j.find("planes") != j.end()) { - return j["planes"]; + if (fileData.find("planes") != fileData.end()) { + return fileData["planes"]; } else { std::array nCells = { - j["numberOfCells"][0], - j["numberOfCells"][1], - j["numberOfCells"][2] + fileData["numberOfCells"][0], + fileData["numberOfCells"][1], + fileData["numberOfCells"][2] }; std::array min, max; - min = j["boundingBox"][0]; - max = j["boundingBox"][1]; + min = fileData["boundingBox"][0]; + max = fileData["boundingBox"][1]; return { utils::GridTools::linspace(min[0], max[0], nCells[0] + 1), @@ -44,48 +47,90 @@ Grid parseGridFromJSON(const nlohmann::json &j) } } -Mesh readMesh(const std::string &fn) +std::vector readObjectsFromJSON(const nlohmann::json& fileData) { - nlohmann::json j; + std::vector objects; - { - std::ifstream i(fn); - i >> j; + if (fileData.contains("objects")) { + for (const auto& obj : fileData["objects"]) { + ObjectDefinition objDef; + objDef.filename = obj["filename"].get(); + objDef.group = obj.value("group", std::filesystem::path(objDef.filename).stem().string()); + if (obj.contains("volume")){ + objDef.isVolume = obj["volume"]; + } + if (obj.contains("mesher")) { + objDef.mesherOverride = obj["mesher"]; + } + objects.push_back(objDef); + } + } else if (fileData.contains("object")) { + ObjectDefinition objDef; + objDef.filename = fileData["object"]["filename"].get(); + objDef.group = std::filesystem::path(objDef.filename).stem().string(); + if (fileData["object"].contains("volume")){ + objDef.isVolume = fileData["object"]["volume"]; + } + if (fileData.contains("mesher")) { + objDef.mesherOverride = fileData["mesher"]; + } + objects.push_back(objDef); + } else { + throw std::runtime_error("No objects defined in input file"); } - std::filesystem::path caseFolder = std::filesystem::path(fn).parent_path(); - std::filesystem::path objPathFromInput = j["object"]["filename"]; - std::filesystem::path meshObjectPath = caseFolder / objPathFromInput; + return objects; +} + +Mesh readMesh(const nlohmann::json& fileData, const std::filesystem::path& folderPath, const ObjectDefinition& objDef) +{ + std::filesystem::path meshObjectPath = folderPath / objDef.filename; std::cout << "-- Reading mesh groups from: " << meshObjectPath; Mesh res = vtkIO::readInputMesh(meshObjectPath); std::cout << "....... [OK]" << std::endl; std::cout << "-- Reading grid from input file"; - res.grid = parseGridFromJSON(j["grid"]); + res.grid = parseGridFromJSON(fileData["grid"]); std::cout << "....... [OK]" << std::endl; + if (res.groups.empty()) { + res.groups.push_back(Group{objDef.group, {}}); + } else { + auto auxResult = utils::meshTools::extractGroupsByName(res, {objDef.group}); + if (auxResult.countElems() != 0) { + return auxResult; + } else { + res.groups[0].name = objDef.group; + } + } + return res; } -std::string readMesherType(const std::string &fn) -{ - nlohmann::json j; +std::string readMesherType(const nlohmann::json& fileData, const std::optional& override) { - std::ifstream i(fn); - i >> j; + nlohmann::json mesherConfig; + + if (override.has_value()) { + mesherConfig = *override; + } else if (fileData.contains("mesher")) { + mesherConfig = fileData["mesher"]; + } else { + return meshlib::app::staircase_mesher; } - if (j["mesher"].contains("type")) { - return j["mesher"]["type"]; + + if (mesherConfig.contains("type")) { + return mesherConfig["type"]; } else { return meshlib::app::staircase_mesher; } } -std::string readExtension(const std::string &fn) +std::string readExtension(const nlohmann::json& fileData, const std::optional& override) { - auto mesherType = readMesherType(fn); + auto mesherType = readMesherType(fileData, override); if (mesherType == meshlib::app::staircase_mesher) { return "str"; } else if (mesherType == meshlib::app::conformal_mesher) { @@ -95,27 +140,79 @@ std::string readExtension(const std::string &fn) } } -meshlib::meshers::ConformalMesherOptions readConformalMesherOptions(const std::string &fn) -{ - nlohmann::json j; - { - std::ifstream i(fn); - i >> j; +meshlib::meshers::StaircaseMesherOptions readStaircaseMesherOptions(const nlohmann::json &fileData, bool isVolume, const std::optional& override) +{ + nlohmann::json mesherConfig; + if (override.has_value()) { + mesherConfig = *override; + } else if (fileData.contains("mesher")) { + mesherConfig = fileData["mesher"]; + } + + meshlib::meshers::StaircaseMesherOptions res; + 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, bool isVolume, const std::optional& override) +{ + nlohmann::json mesherConfig; + if (override.has_value()) { + mesherConfig = *override; + } else if (fileData.contains("mesher")) { + mesherConfig = fileData["mesher"]; + } + meshlib::meshers::ConformalMesherOptions res; - if (j["mesher"].contains("options")) { - res.snapperOptions.edgePoints = j["mesher"]["options"]["edgePoints"]; - res.snapperOptions.forbiddenLength = j["mesher"]["options"]["forbiddenLength"]; + if (isVolume){ + res.volumeGroups.insert(0); + } + if (mesherConfig.contains("options")) { + res.snapperOptions.edgePoints = mesherConfig["options"]["edgePoints"]; + res.snapperOptions.forbiddenLength = mesherConfig["options"]["forbiddenLength"]; } return res; } -std::unique_ptr buildMesher(const Mesh &in, const std::string &fn) + +bool readExportGridOption(const nlohmann::json& fileData, const std::optional& override) +{ + nlohmann::json mesherConfig; + if (override.has_value()) { + mesherConfig = *override; + } else if (fileData.contains("mesher")) { + mesherConfig = fileData["mesher"]; + } + + if (mesherConfig.contains("options") && + mesherConfig["options"].contains("exportGrid")) { + return mesherConfig["options"]["exportGrid"]; + } + return true; +} + +std::unique_ptr buildMesher(const Mesh& in, const nlohmann::json & fileData, const ObjectDefinition& objDef) { - auto mesherType = readMesherType(fn); + auto mesherType = readMesherType(fileData, objDef.mesherOverride); + if (mesherType == meshlib::app::staircase_mesher) { - return std::make_unique(meshlib::meshers::StaircaseMesher{in}); + return std::make_unique(meshlib::meshers::StaircaseMesher{ + in, + 4, + readStaircaseMesherOptions(fileData, objDef.isVolume, objDef.mesherOverride) + }); } else if (mesherType == meshlib::app::conformal_mesher) { - return std::make_unique(meshlib::meshers::ConformalMesher{in, readConformalMesherOptions(fn)}); + return std::make_unique(meshlib::meshers::ConformalMesher{in, readConformalMesherOptions(fileData, objDef.isVolume, objDef.mesherOverride)}); } else { throw std::runtime_error("Unsupported mesher type"); } @@ -138,25 +235,47 @@ int launcher(int argc, const char* argv[]) return EXIT_SUCCESS; } - // Input - std::string inputFilename = vm["input"].as(); - std::cout << "-- Input file is: " << inputFilename << std::endl; + std::string inputFileName = vm["input"].as(); + std::cout << "-- Input file is: " << inputFileName << std::endl; - Mesh mesh = readMesh(inputFilename); + nlohmann::json inputFileData; + { + std::ifstream i(inputFileName); + i >> inputFileData; + } + std::vector objects = readObjectsFromJSON(inputFileData); + std::filesystem::path outputFolder = getFolder(inputFileName); + auto basename = getBasename(inputFileName); - // Mesh - auto mesher = buildMesher(mesh, inputFilename); - Mesh resultMesh = mesher->mesh(); + Mesh firstMesh; + bool first = true; - std::filesystem::path outputFolder = getFolder(inputFilename); - auto basename = getBasename(inputFilename); - auto extension = readExtension(inputFilename); - - exportMeshToVTU(outputFolder / (basename + ".tessellator." + extension + ".vtk"), resultMesh); - exportGridToVTU(outputFolder / (basename + ".tessellator.grid.vtk"), resultMesh.grid); + for (const auto& objDef : objects) { + std::cout << "\n-- Processing object: " << objDef.filename << " (group: " << objDef.group << ")" << std::endl; + + Mesh mesh = readMesh(inputFileData, outputFolder, objDef); + + auto mesher = buildMesher(mesh, inputFileData, objDef); + Mesh resultMesh = mesher->mesh(); + + if (first) { + firstMesh = resultMesh; + first = false; + } + + auto extension = readExtension(inputFileData, objDef.mesherOverride); + std::string outputFileName = objDef.group + ".tessellator." + extension + ".vtk"; + exportMeshToVTU(outputFolder / outputFileName, resultMesh); + std::cout << "-- Exported: " << outputFileName << std::endl; + } + + if (!first && readExportGridOption(inputFileData, std::nullopt)) { + exportGridToVTU(outputFolder / (basename + ".tessellator.grid.vtk"), firstMesh.grid); + std::cout << "-- Exported grid: " << basename << ".tessellator.grid.vtk" << std::endl; + } return EXIT_SUCCESS; } -} \ No newline at end of file +} diff --git a/src/app/launcher.h b/src/app/launcher.h index d65d218..48823c3 100644 --- a/src/app/launcher.h +++ b/src/app/launcher.h @@ -1,7 +1,10 @@ #pragma once #include "types/Mesh.h" +#include "meshers/MesherBase.h" #include +#include +#include #include namespace meshlib::app { @@ -9,7 +12,17 @@ namespace meshlib::app { const std::string conformal_mesher ("conformal"); const std::string staircase_mesher ("staircase"); +struct ObjectDefinition { + std::string filename; + std::string group; + bool isVolume = false; + std::optional mesherOverride; +}; + int launcher(int argc, const char* argv[]); -Grid parseGridFromJSON(const nlohmann::json& j); +Grid parseGridFromJSON(const nlohmann::json& fileData); +std::vector readObjectsFromJSON(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 9457e14..71f00b5 100644 --- a/src/app/vtkIO.cpp +++ b/src/app/vtkIO.cpp @@ -1,13 +1,16 @@ #include "vtkIO.h" -#include #include #include #include +#include +#include #include +#include #include #include #include +#include #include #include @@ -43,7 +46,7 @@ vtkSmartPointer readAsVTU(const std::filesystem::path& file } vtkSmartPointer vtu; - std::string extension = vtksys::SystemTools::GetFilenameLastExtension(fn); + std::string extension = fn.substr(fn.find_last_of(".")).empty() ? "" : fn.substr(fn.find_last_of(".")); std::transform(extension.begin(), extension.end(), extension.begin(), ::tolower); @@ -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( @@ -163,12 +189,28 @@ vtkSmartPointer toVTKGroupsArray(const Mesh& mesh) return groupsDataArray; } +vtkSmartPointer toVTKGroupNamesArray(const Mesh& mesh) +{ + vtkNew groupNamesArray; + groupNamesArray->SetName("groupNames"); + groupNamesArray->SetNumberOfComponents(1); + + for (const auto& group : mesh.groups) { + for (std::size_t e = 0; e < group.elements.size(); e++) { + groupNamesArray->InsertNextValue(group.name.c_str()); + } + } + + return groupNamesArray; +} + vtkSmartPointer elementsToVTU(const Mesh& mesh) { vtkNew vtu; vtu->SetPoints(toVTKPoints(mesh.coordinates)); vtu->GetCellData()->AddArray(toVTKGroupsArray(mesh)); + vtu->GetCellData()->AddArray(toVTKGroupNamesArray(mesh)); std::vector cellTypes; cellTypes.reserve(mesh.countElems()); @@ -189,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 d481f43..b3dc8a3 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -2,11 +2,15 @@ message(STATUS "Creating build system for tessellator-core") add_library(tessellator-core "Collapser.cpp" + "Compressor.cpp" "Slicer.cpp" "Snapper.cpp" "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/Compressor.cpp b/src/core/Compressor.cpp new file mode 100644 index 0000000..b8c84ca --- /dev/null +++ b/src/core/Compressor.cpp @@ -0,0 +1,765 @@ +#include "Compressor.h" + +#include +#include + +#include "types/Mesh.h" +#include "utils/GridTools.h" + +namespace meshlib::core { + +using meshlib::Sign; +using meshlib::SignedAxis; +using meshlib::PlanePoint; +using meshlib::PlaneLinel; +using meshlib::PlaneSurfel; +using meshlib::PlaneSurface; +using meshlib::Contour; +using meshlib::CrossLine; + +std::size_t Compressor::compressSurfacesInMesh(Mesh& mesh) { + std::size_t totalOriginal = 0; + std::size_t totalCompressed = 0; + + for (Group& group : mesh.groups) { + std::vector surfaces; + + for (const Element& elem : group.elements) { + if (elem.type == Element::Type::Surface) { + surfaces.push_back(elem); + } + } + + if (surfaces.empty()) { + continue; + } + + totalOriginal += surfaces.size(); + std::vector compressedSurfaces = compressSurfaces_(mesh.coordinates, surfaces); + totalCompressed += compressedSurfaces.size(); + + // Build new elements vector with compressed surfaces + std::vector newElements; + ElementId surfaceIdx = 0; + for (ElementId e = 0; e < group.elements.size(); e++) { + if (group.elements[e].type == Element::Type::Surface) { + if (surfaceIdx < compressedSurfaces.size()) { + newElements.push_back(compressedSurfaces[surfaceIdx]); + surfaceIdx++; + } + } else { + newElements.push_back(group.elements[e]); + } + } + group.elements = std::move(newElements); + } + + return totalOriginal - totalCompressed; +} + +std::size_t Compressor::compressLinesInMesh(Mesh& mesh, const std::vector & dimensionPolicy) { + std::size_t totalOriginal = 0; + std::size_t totalCompressed = 0; + + std::vector compress; + if (dimensionPolicy.size() == 0){ + compress = std::vector(mesh.groups.size(), true); + } + else{ + compress.reserve(mesh.groups.size()); + for (auto policy: dimensionPolicy){ + compress.emplace_back(policy == Element::Type::Volume || policy == Element::Type::Surface); + } + } + + for (GroupId g = 0; g < mesh.groups.size(); ++g){ + if (!compress[g]){ + continue; + } + + std::vector lines; + + for (const Element& elem : mesh.groups[g].elements) { + if (elem.type == Element::Type::Line) { + lines.push_back(elem); + } + } + + if (lines.empty()) { + continue; + } + + totalOriginal += lines.size(); + std::vector compressedLines = compressLines_(mesh.coordinates, lines); + totalCompressed += compressedLines.size(); + + // Build new elements vector with compressed lines + std::vector newElements; + ElementId lineIdx = 0; + for (ElementId e = 0; e < mesh.groups[g].elements.size(); e++) { + if (mesh.groups[g].elements[e].type == Element::Type::Line) { + if (lineIdx < compressedLines.size()) { + newElements.push_back(compressedLines[lineIdx]); + lineIdx++; + } + } else { + newElements.push_back(mesh.groups[g].elements[e]); + } + } + mesh.groups[g].elements = std::move(newElements); + } + + return totalOriginal - totalCompressed; +} + +std::vector Compressor::compressLines_( + const std::vector& coords, + const std::vector& lines) { + std::vector result; + std::map, + std::vector> signDirLines; + for (std::size_t l = 0; l < lines.size(); l++) { + std::array auxCells; + auxCells[0] = utils::GridTools::toCell(coords[lines[l].vertices[0]]); + auxCells[1] = utils::GridTools::toCell(coords[lines[l].vertices[1]]); + GridLine gridLine; + Sign sign = 1; + Axis direction = 0; + for (Axis d = X; X <= Z; d++) { + if (auxCells[0](d) != auxCells[1](d)) { + if (auxCells[0](d) > auxCells[1](d)) { + sign = -1; + } + direction = d; + Axis d1 = (d + 1) % 3; + Axis d2 = (d + 2) % 3; + gridLine[0] = auxCells[0](d1); + gridLine[1] = auxCells[0](d2); + break; + } + } + signDirLines[std::make_pair(gridLine, + std::make_pair(sign, direction))].push_back(l); + } + for (std::map, + std::vector>::const_iterator + it = signDirLines.begin(); it != signDirLines.end(); ++it) { + std::vector auxElems; + for (std::size_t i = 0; i < it->second.size(); i++) { + auxElems.push_back(lines[it->second[i]]); + } + std::vector compressedLines = + compressDirSignLines_(coords, + it->first.second, + auxElems); + result.insert(result.end(), compressedLines.begin(), compressedLines.end()); + } + return result; +} + +std::vector Compressor::compressDirSignLines_( + const std::vector& coords, + const SignedAxis& signedDir, + const std::vector& lines) { + std::vector result; + std::map> relativeLines; + std::map> lineRelatives; + for (std::size_t l = 0; l < lines.size(); l++) { + for (RelativeId vertex : lines[l].vertices) { + relativeLines[vertex].insert(l); + lineRelatives[l].insert(vertex); + } + } + std::set visitedLineIds; + for (std::map>::const_iterator + itExt = lineRelatives.begin(); itExt != lineRelatives.end(); ++itExt) { + if (visitedLineIds.count(itExt->first) == 0) { + RelativeId minCell = lines[itExt->first].vertices[0]; + RelativeId maxCell = lines[itExt->first].vertices[1]; + std::queue linesToVisit; + linesToVisit.push(itExt->first); + visitedLineIds.insert(itExt->first); + while (!linesToVisit.empty()) { + ElementId elem = linesToVisit.front(); + linesToVisit.pop(); + for (std::size_t i = 0; i < 2; i++) { + if (coords[minCell] > coords[lines[elem].vertices[i]]) { + minCell = lines[elem].vertices[i]; + } + if (coords[maxCell] < coords[lines[elem].vertices[i]]) { + maxCell = lines[elem].vertices[i]; + } + } + for (std::set::const_iterator + itCell = lineRelatives[elem].begin(); + itCell != lineRelatives[elem].end(); ++itCell) { + for (std::set::const_iterator + itLine = relativeLines[*itCell].begin(); + itLine != relativeLines[*itCell].end(); ++itLine) { + if (visitedLineIds.count(*itLine) == 0) { + linesToVisit.push(*itLine); + visitedLineIds.insert(*itLine); + } + } + } + } + Element newElem; + newElem.type = Element::Type::Line; + newElem.vertices.push_back(minCell); + newElem.vertices.push_back(maxCell); + if (signedDir.first < 0) { + std::swap(newElem.vertices[0], newElem.vertices[1]); + } + result.push_back(newElem); + } + } + return result; +} + +std::vector Compressor::compressSurfaces_( + std::vector& coords, + const std::vector& surfaces) { + + std::vector result; + std::map, std::vector> signDirSurfs; + + for (std::size_t s = 0; s < surfaces.size(); s++) { + if (surfaces[s].vertices.size() != 4) { + result.push_back(surfaces[s]); + continue; + } + std::array auxCells; + auxCells[0] = utils::GridTools::toCell(coords[surfaces[s].vertices[0]]); + auxCells[1] = utils::GridTools::toCell(coords[surfaces[s].vertices[1]]); + auxCells[2] = utils::GridTools::toCell(coords[surfaces[s].vertices[2]]); + CellDir gridSurface; + Sign sign = 1; + Axis direction = 0; + for (Axis d = X; d <= Z; d++) { + if (auxCells[0](d) == auxCells[2](d)) { + Cell normal = (auxCells[1] - auxCells[0]) ^ + (auxCells[2] - auxCells[0]); + if (normal(d) >= 0) { + sign = 1; + } else { + sign = -1; + } + direction = d; + gridSurface = auxCells[0](d); + break; + } + } + signDirSurfs[std::make_pair(gridSurface, + std::make_pair(sign, direction))].push_back(s); + } + for (std::map, + std::vector>::const_iterator + it = signDirSurfs.begin(); it != signDirSurfs.end(); ++it) { + std::vector auxElems; + for (std::size_t i = 0; i < it->second.size(); i++) { + auxElems.push_back(surfaces[it->second[i]]); + } + std::vector compressedSurfaces = + compressSurfacesWithSameNormal_(coords, it->first.second, auxElems); + result.insert(result.end(), compressedSurfaces.begin(), compressedSurfaces.end()); + } + return result; +} + +std::vector Compressor::compressSurfacesWithSameNormal_( + std::vector& coords, + const SignedAxis& signedDir, + const std::vector& surfaces) { + std::vector result; + std::map> edgeSurfaces; + std::map> surfaceEdges; + for (ElementId s = 0; s < surfaces.size(); s++) { + for (std::size_t i = 0; i < 4; i++) { + std::size_t j = (i + 1) % 4; + LinIds edge; + edge[0] = surfaces[s].vertices[i]; + edge[1] = surfaces[s].vertices[j]; + std::sort(edge.begin(), edge.end()); + edgeSurfaces[edge].insert(s); + surfaceEdges[s].insert(edge); + } + } + std::set visitedSurfaceIds; + for (std::map>::const_iterator + itExt = surfaceEdges.begin(); itExt != surfaceEdges.end(); ++itExt) { + if (visitedSurfaceIds.count(itExt->first) == 0) { + std::set connectedSurfaceIds; + std::queue surfacesToVisit; + surfacesToVisit.push(itExt->first); + visitedSurfaceIds.insert(itExt->first); + while (!surfacesToVisit.empty()) { + ElementId e = surfacesToVisit.front(); + surfacesToVisit.pop(); + connectedSurfaceIds.insert(e); + for (std::set::const_iterator + itLine = surfaceEdges[e].begin(); + itLine != surfaceEdges[e].end(); ++itLine) { + for (std::set::const_iterator + itSurf = edgeSurfaces[*itLine].begin(); + itSurf != edgeSurfaces[*itLine].end(); ++itSurf) { + if (visitedSurfaceIds.count(*itSurf) == 0) { + surfacesToVisit.push(*itSurf); + visitedSurfaceIds.insert(*itSurf); + } + } + } + } + std::vector connectedSurfaces; + for (std::set::const_iterator + it = connectedSurfaceIds.begin(); it != connectedSurfaceIds.end(); ++it) { + connectedSurfaces.push_back(surfaces[*it]); + } + std::vector compressedSurfaces = compressConnectedSurfaces_(coords, signedDir, connectedSurfaces); + result.insert(result.end(), compressedSurfaces.begin(), compressedSurfaces.end()); + } + } + return result; +} + +std::vector Compressor::compressConnectedSurfaces_( + std::vector& coords, + const SignedAxis& signDir, + const std::vector& surfs) { + std::vector result; + Axis d = signDir.second; + Axis d1 = (d + 1) % 3; + Axis d2 = (d + 2) % 3; + CellDir plane = utils::GridTools::toCell(coords[surfs[0].vertices[0]])(d); + std::set surfels; + for (std::size_t s = 0; s < surfs.size(); s++) { + std::pair fullSurfacePoints; + fullSurfacePoints.first[0] = + utils::GridTools::toCell(coords[surfs[s].vertices[0]])(d1); + fullSurfacePoints.first[1] = + utils::GridTools::toCell(coords[surfs[s].vertices[0]])(d2); + fullSurfacePoints.second[0] = + utils::GridTools::toCell(coords[surfs[s].vertices[2]])(d1); + fullSurfacePoints.second[1] = + utils::GridTools::toCell(coords[surfs[s].vertices[2]])(d2); + CellDir i0 = std::min(fullSurfacePoints.first[0], fullSurfacePoints.second[0]); + CellDir i1 = std::max(fullSurfacePoints.first[0], fullSurfacePoints.second[0]); + CellDir j0 = std::min(fullSurfacePoints.first[1], fullSurfacePoints.second[1]); + CellDir j1 = std::max(fullSurfacePoints.first[1], fullSurfacePoints.second[1]); + for (CellDir i = i0; i < i1; i++) { + for (CellDir j = j0; j < j1; j++) { + PlaneSurfel surfel = {{i, j}}; + surfels.insert(surfel); + } + } + } + std::vector maximalRectangles = compressSurfelsIntoMaximalRectangles_(surfels); + CoordinateMap coordMap; + for (std::size_t s = 0; s < surfs.size(); s++) { + for (std::size_t i = 0; i < 4; i++) { + RelativeId coordId = surfs[s].vertices[i]; + Relative coord = coords[coordId]; + coordMap[coord] = coordId; + } + } + for (std::size_t e = 0; e < maximalRectangles.size(); e++) { + std::array corners; + corners[0](d) = corners[2](d) = plane; + corners[0](d1) = maximalRectangles[e].first[0]; + corners[0](d2) = maximalRectangles[e].first[1]; + corners[2](d1) = maximalRectangles[e].second[0]; + corners[2](d2) = maximalRectangles[e].second[1]; + corners[1] = corners[3] = corners[0]; + corners[1](d1) = corners[2](d1); + corners[3](d2) = corners[2](d2); + if (signDir.first < 0) { + std::swap(corners[1], corners[3]); + } + Element newSurface; + newSurface.type = Element::Type::Surface; + for (std::size_t i = 0; i < 4; i++) { + Relative rel = utils::GridTools::toRelative(corners[i]); + if (coordMap.count(rel) == 0) { + coordMap[rel] = coords.size(); + coords.push_back(rel); + } + newSurface.vertices.push_back(coordMap[rel]); + } + result.push_back(newSurface); + } + return result; +} + +std::vector Compressor::compressSurfelsIntoMaximalRectangles_( + const std::set& surfels) { + std::vector result; + const std::vector& contours = getContours_(surfels); + std::array, 2> crossingLines = getCrossingLines_(surfels, contours); + crossingLines = getMaxCompactedLines_(crossingLines); + std::set linels; + for (std::size_t c = 0; c < contours.size(); c++) { + for (std::size_t i = 0; i < contours[c].size(); i++) { + std::size_t j = (i + 1) % contours[c].size(); + std::set aux = getLinelsBetween_(contours[c][i], contours[c][j]); + linels.insert(aux.begin(), aux.end()); + } + } + for (Axis d = 0; d < 2; d++) { + Axis d1 = (d + 1) % 2; + for (std::size_t i = 0; i < crossingLines[d].size(); i++) { + PlanePoint ini, end; + ini[d1] = end[d1] = crossingLines[d][i].first; + ini[d] = crossingLines[d][i].second.first; + end[d] = crossingLines[d][i].second.second; + std::set aux = getLinelsBetween_(ini, end); + linels.insert(aux.begin(), aux.end()); + } + } + addConcaveLinels_(surfels, contours, linels); + std::map> edgeSurfels; + std::map> surfelEdges; + for (std::set::const_iterator + it = surfels.begin(); it != surfels.end(); ++it) { + surfelEdges.insert(std::make_pair(*it, std::set())); + for (Axis d = 0; d < 2; d++) { + for (CellDir diff = -1; diff <= 1; diff += 2) { + PlaneLinel linel = getSurfaceEdge_(*it, diff, d); + if (linels.count(linel) == 0) { + surfelEdges[*it].insert(linel); + edgeSurfels[linel].insert(*it); + } + } + } + } + std::set visitedSurfels; + for (std::map>::const_iterator + itSurfExt = surfelEdges.begin(); + itSurfExt != surfelEdges.end(); ++itSurfExt) { + if (visitedSurfels.count(itSurfExt->first) == 0) { + std::queue surfelsToVisit; + surfelsToVisit.push(itSurfExt->first); + visitedSurfels.insert(itSurfExt->first); + PlanePoint minPoint = itSurfExt->first; + PlanePoint maxPoint = itSurfExt->first; + while (!surfelsToVisit.empty()) { + PlaneSurfel surfel = surfelsToVisit.front(); + surfelsToVisit.pop(); + if (surfel < minPoint) { + minPoint = surfel; + } + if (surfel > maxPoint) { + maxPoint = surfel; + } + for (std::set::const_iterator + itLin = surfelEdges[surfel].begin(); + itLin != surfelEdges[surfel].end(); ++itLin) { + for (std::set::const_iterator + itSurfInt = edgeSurfels[*itLin].begin(); + itSurfInt != edgeSurfels[*itLin].end(); ++itSurfInt) { + if (visitedSurfels.count(*itSurfInt) == 0) { + surfelsToVisit.push(*itSurfInt); + visitedSurfels.insert(*itSurfInt); + } + } + } + } + maxPoint[0]++; + maxPoint[1]++; + result.push_back(std::make_pair(minPoint, maxPoint)); + } + } + return result; +} + +std::vector Compressor::getContours_(const std::set& surfels) { + std::vector result; + if (surfels.empty()) { + return result; + } + std::set visitedEdges; + result.push_back(getContourFromStartingEdge_( + getSurfaceEdge_(*surfels.begin(), -1, 0), + surfels, + visitedEdges + )); + for (std::set::const_iterator + it = surfels.begin(); it != surfels.end(); ++it) { + for (Axis d = 0; d < 2; d++) { + for (CellDir diff = -1; diff <= 1; diff += 2) { + PlaneSurfel adjSurf; + PlaneLinel adjEdge; + adjSurf = *it; + adjSurf[d] += diff; + adjEdge = getSurfaceEdge_(*it, diff, d); + if ((surfels.find(adjSurf) == surfels.end()) && + (visitedEdges.find(adjEdge) == visitedEdges.end())) { + result.push_back(getContourFromStartingEdge_(adjEdge, surfels, visitedEdges)); + } + } + } + } + return result; +} + +Contour Compressor::getContourFromStartingEdge_( + const PlaneLinel& from, + const std::set& surfs, + std::set& visitedEdges) { + Contour result; + std::queue q; + if (visitedEdges.find(from) != visitedEdges.end()) { + return result; + } + std::vector lines; + q.push(from); + visitedEdges.insert(from); + lines.push_back(from); + while (!q.empty()) { + PlaneLinel edge = q.front(); + q.pop(); + PlanePoint pos = edge.first; + Axis d0 = edge.second; + Axis d1 = (d0 + 1) % 2; + PlaneSurfel surf = pos; + if (surfs.find(surf) == surfs.end()) { + surf[d1]--; + } + for (CellDir diff = -1; diff <= 1; diff += 2) { + PlaneSurfel adjSurf1 = surf; + adjSurf1[d0] += diff; + if (surfs.find(adjSurf1) == surfs.end()) { + PlaneLinel adjEdge = getSurfaceEdge_(surf, diff, d0); + if (visitedEdges.find(adjEdge) == visitedEdges.end()) { + q.push(adjEdge); + visitedEdges.insert(adjEdge); + lines.push_back(adjEdge); + break; + } + continue; + } + if (surf == pos) { + adjSurf1[d1]--; + } else { + adjSurf1[d1]++; + } + if (surfs.find(adjSurf1) == surfs.end()) { + PlaneLinel adjEdge = edge; + adjEdge.first[d0] += diff; + if (visitedEdges.find(adjEdge) == visitedEdges.end()) { + q.push(adjEdge); + visitedEdges.insert(adjEdge); + lines.push_back(adjEdge); + break; + } + continue; + } else { + PlaneLinel adjEdge = getSurfaceEdge_(adjSurf1, -diff, d0); + if (visitedEdges.find(adjEdge) == visitedEdges.end()) { + q.push(adjEdge); + visitedEdges.insert(adjEdge); + lines.push_back(adjEdge); + break; + } + continue; + } + } + } + for (std::vector::const_iterator + it = lines.begin(); it != lines.end(); ++it) { + std::vector::const_iterator itPlus = std::next(it); + if (itPlus == lines.end()) { + itPlus = lines.begin(); + } + if (it->second == itPlus->second) { + continue; + } + std::array extremes = {it->first, it->first}; + extremes[1][it->second]++; + std::array extremesP = {itPlus->first, itPlus->first}; + extremesP[1][itPlus->second]++; + for (std::size_t p = 0; p < 4; p++) { + if (extremes[p / 2] == extremesP[p % 2]) { + result.push_back(extremes[p / 2]); + break; + } + } + } + return result; +} + +PlaneLinel Compressor::getSurfaceEdge_(const PlaneSurfel& surfel, const CellDir& diff, const Axis& dir) { + PlaneLinel result; + result.first = surfel; + result.second = (dir + 1) % 2; + if (diff > 0) { + result.first[dir]++; + } + return result; +} + +std::array, 2> + Compressor::getCrossingLines_( + const std::set& surfs, + const std::vector& conts) { + std::array, 2> result; + for (Axis d = 0; d < 2; d++) { + Axis d1 = (d + 1) % 2; + std::map> cross; + for (std::vector::const_iterator + it1 = conts.begin(); it1 != conts.end(); ++it1) { + for (std::vector::const_iterator + it2 = it1->begin(); it2 != it1->end(); ++it2) { + cross[(*it2)[d1]].insert((*it2)[d]); + } + } + for (std::map>::const_iterator + itMap = cross.begin(); itMap != cross.end(); ++itMap) { + for (std::set::const_iterator + itSet = itMap->second.begin(); + itSet != itMap->second.end(); ++itSet) { + std::set::const_iterator itSetPlus = std::next(itSet); + if (itSetPlus == itMap->second.end()) { + break; + } + bool valid = true; + PlanePoint edge; + edge[d1] = itMap->first; + for (CellDir i = *itSet; i < *itSetPlus; i++) { + edge[d] = i; + PlaneSurfel adjSurf1, adjSurf2; + adjSurf1 = adjSurf2 = edge; + adjSurf1[d1]--; + if ((surfs.find(adjSurf1) == surfs.end()) || + (surfs.find(adjSurf2) == surfs.end())) { + valid = false; + break; + } + } + if (valid) { + result[d].push_back( + std::make_pair(itMap->first, + std::make_pair(*itSet, *itSetPlus))); + } + } + } + } + return result; +} + +std::array, 2> + Compressor::getMaxCompactedLines_( + const std::array, 2>& cross) { + std::array, 2> result; + if (cross[0].size() > cross[1].size()) { + result[0] = cross[0]; + } else { + result[1] = cross[1]; + } + return result; +} + +std::set Compressor::getLinelsBetween_( + const PlanePoint& ini, + const PlanePoint& end) { + std::set result; + for (Axis d = 0; d < 2; d++) { + Axis d1 = (d + 1) % 2; + if (ini[d] == end[d]) { + PlaneLinel linel; + linel.first[d] = ini[d]; + linel.second = d1; + for (CellDir + k = std::min(ini[d1], end[d1]); + k < std::max(ini[d1], end[d1]); k++) { + linel.first[d1] = k; + result.insert(linel); + } + } + } + return result; +} + +void Compressor::addConcaveLinels_(const std::set& surfs, + const std::vector& conts, + std::set& lines) { + std::set concavePoints; + for (std::size_t c = 0; c < conts.size(); c++) { + for (std::size_t i = 0; i < conts[c].size(); i++) { + PlanePoint point = conts[c][i]; + std::size_t numSurfAdj = 0; + std::size_t numLineAdj = 0; + for (CellDir diffx = -1; diffx < 1; diffx++) { + for (CellDir diffy = -1; diffy < 1; diffy++) { + PlaneSurfel surfel = point; + surfel[0] += diffx; + surfel[1] += diffy; + if (surfs.count(surfel) != 0) { + numSurfAdj++; + } + } + } + for (Axis d = 0; d < 2; d++) { + for (CellDir diff = -1; diff < 1; diff++) { + PlaneLinel linel = std::make_pair(point, d); + linel.first[d] += diff; + if (lines.count(linel) != 0) { + numLineAdj++; + } + } + } + if ((numSurfAdj > 2) && (numLineAdj < 3)) { + concavePoints.insert(point); + } + } + } + for (std::set::const_iterator + it = concavePoints.begin(); it != concavePoints.end(); ++it) { + std::size_t numLineAdj = 0; + for (Axis d = 0; d < 2; d++) { + for (CellDir diff = -1; diff < 1; diff++) { + PlaneLinel linel = std::make_pair(*it, d); + linel.first[d] += diff; + if (lines.count(linel) != 0) { + numLineAdj++; + } + } + } + if (numLineAdj > 2) { + continue; + } + for (Axis d = 0; d < 2; d++) { + Axis d1 = (d + 1) % 2; + bool found = false; + for (CellDir diff = -1; diff < 1; diff++) { + PlaneLinel linel = std::make_pair(*it, d); + linel.first[d] += diff; + if (lines.count(linel) == 0) { + lines.insert(linel); + while (true) { + PlaneLinel aux1, aux2; + if (diff < 0) { + aux1 = aux2 = std::make_pair(linel.first, d1); + aux1.first[d1]--; + linel.first[d]--; + } else { + linel.first[d]++; + aux1 = aux2 = std::make_pair(linel.first, d1); + aux1.first[d1]--; + } + if ((lines.count(aux1) != 0) || + (lines.count(aux2) != 0)) { + break; + } + lines.insert(linel); + } + found = true; + break; + } + } + if (found) { + break; + } + } + } +} + +} diff --git a/src/core/Compressor.h b/src/core/Compressor.h new file mode 100644 index 0000000..57650be --- /dev/null +++ b/src/core/Compressor.h @@ -0,0 +1,90 @@ +#pragma once + +#include "types/Mesh.h" +#include "utils/Types.h" + +#include +#include + +namespace meshlib::core { + +class Compressor { +public: + // Compress only quad surfaces (4 vertices) that are coplanar and adjacent + // Returns number of surfaces merged (original_count - compressed_count) + static std::size_t compressSurfacesInMesh(Mesh& mesh); + + // Compress collinear line segments that are adjacent + // Returns number of lines merged (original_count - compressed_count) + static std::size_t compressLinesInMesh(Mesh& mesh, const std::vector & dimensionPolicy = {}); + +private: + // Group surfaces by (grid_plane, sign, axis) and compress each group + static std::vector compressSurfaces_( + std::vector& coords, + const std::vector& surfs); + + // Compress surfaces with same normal direction and sign + static std::vector compressSurfacesWithSameNormal_( + std::vector& coords, + const std::pair& signDir, + const std::vector& surfs); + + // Compress connected coplanar surfaces using contour detection + static std::vector compressConnectedSurfaces_( + std::vector& coords, + const std::pair& signDir, + const std::vector& surfs); + + // Group lines by (grid_line, sign, axis) and compress each group + static std::vector compressLines_( + const std::vector& coords, + const std::vector& lines); + + // Compress lines with same direction and sign + static std::vector compressDirSignLines_( + const std::vector& coords, + const std::pair& signDir, + const std::vector& lines); + + // Merge adjacent surfels into maximal rectangles + static std::vector compressSurfelsIntoMaximalRectangles_( + const std::set& surfs); + + // Detect boundary contours of surfel set + static std::vector getContours_(const std::set& surfs); + + // Trace a single contour from starting edge + static Contour getContourFromStartingEdge_( + const PlaneLinel& from, + const std::set& surfs, + std::set& visited); + + // Get edge of a surfel in given direction + static PlaneLinel getSurfaceEdge_( + const PlaneSurfel& surf, + const CellDir& diff, + const Axis& dir); + + // Find lines crossing between contours + static std::array, 2> getCrossingLines_( + const std::set& surfs, + const std::vector& contours); + + // Select the set of crossing lines with maximum count + static std::array, 2> getMaxCompactedLines_( + const std::array, 2>& cross); + + // Get linels between two points + static std::set getLinelsBetween_( + const PlanePoint& ini, + const PlanePoint& end); + + // Add linels at concave corners + static void addConcaveLinels_( + const std::set& surfs, + const std::vector& contours, + std::set& lines); +}; + +} 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/Smoother.cpp b/src/core/Smoother.cpp index 9a15452..9dc38ea 100644 --- a/src/core/Smoother.cpp +++ b/src/core/Smoother.cpp @@ -8,6 +8,8 @@ #include #include +#include +#include #ifdef TESSELLATOR_EXECUTION_POLICIES #include @@ -19,6 +21,63 @@ namespace core { using namespace utils; using namespace meshTools; +namespace { + +using GridEntity = std::pair; + +struct OwnedCoordinates { + std::map edges; + std::map faces; + std::map interiors; +}; + +Cell canonicalCell(Cell cell, const SmootherTools& tools) +{ + for (Axis axis = X; axis <= Z; ++axis) { + cell[axis] = std::max(0, std::min(cell[axis], tools.numCellsDir(axis) - 1)); + } + return cell; +} + +OwnedCoordinates buildOwnedCoordinates( + const Coordinates& coordinates, + const SmootherTools& tools) +{ + OwnedCoordinates owned; + for (CoordinateId id = 0; id < coordinates.size(); ++id) { + const auto& coordinate = coordinates[id]; + if (tools.isRelativeInCellCorner(coordinate)) { + continue; // Grid corners are fixed throughout smoothing. + } + + Cell cell = canonicalCell(tools.toCell(coordinate), tools); + const auto edge = tools.getCellEdgeAxis(coordinate); + if (edge.first) { + owned.edges[{cell, edge.second}].insert(id); + continue; + } + const auto face = tools.getCellFaceAxis(coordinate); + if (face.first) { + owned.faces[{cell, face.second}].insert(id); + continue; + } + owned.interiors[cell].insert(id); + } + return owned; +} + +SmootherTools::IncidentElements buildIncidentElements(const Elements& elements) +{ + SmootherTools::IncidentElements incident; + for (const auto& element : elements) { + for (const auto id : element.vertices) { + incident[id].push_back(&element); + } + } + return incident; +} + +} Smoother::Smoother(const Mesh& mesh, const SmootherOptions& opts) : sT_(SmootherTools(mesh.grid)), @@ -35,6 +94,9 @@ Smoother::Smoother(const Mesh& mesh, const SmootherOptions& opts) : auto const singularIds = sT_.buildSingularIds(g.elements, mesh_.coordinates, opts_.featureDetectionAngle); + // Boundary remeshing is preprocessing. Keep its smooth sets separate + // from the ownership phases below: an element can touch several grid + // entities, while a coordinate has exactly one canonical owner. std::vector patchs; for (auto const& cell : sT_.buildCellElemMap(g.elements, mesh_.coordinates)) { for (auto const& p : @@ -46,34 +108,45 @@ Smoother::Smoother(const Mesh& mesh, const SmootherOptions& opts) : std::for_each(patchs.begin(), patchs.end(), [&](auto& p) { sT_.remeshBoundary(g.elements, res.coordinates, mesh_.coordinates, p); }); - - std::for_each(patchs.begin(), patchs.end(), [&](auto& p) { - sT_.collapsePointsOnCellEdges(res.coordinates, p, singularIds, opts_.contourAlignmentAngle); - }); - std::for_each( -#ifdef TESSELLATOR_EXECUTION_POLICIES - std::execution::par, -#endif - patchs.begin(), patchs.end(), [&](auto& p) { - sT_.collapsePointsOnCellFaces(res.coordinates, p, singularIds); - }); + const auto incidentElements = buildIncidentElements(g.elements); + const auto owned = buildOwnedCoordinates(res.coordinates, sT_); + IdSet edgeIds, faceIds, interiorIds; + for (const auto& entity : owned.edges) { + edgeIds.insert(entity.second.begin(), entity.second.end()); + } + for (const auto& entity : owned.faces) { + faceIds.insert(entity.second.begin(), entity.second.end()); + } + for (const auto& entity : owned.interiors) { + interiorIds.insert(entity.second.begin(), entity.second.end()); + } - std::for_each( -#ifdef TESSELLATOR_EXECUTION_POLICIES - std::execution::par, -#endif - patchs.begin(), patchs.end(), [&](auto& p) { - sT_.collapsePointsOnFeatureEdges(res.coordinates, p, singularIds); - }); + // Later phases may read an earlier boundary but never move it. + // This is intentionally sequential: entity ownership makes a future + // conflict graph possible, but entities sharing an element still need + // serialization until that scheduler is introduced. + for (const auto& patch : patchs) { + sT_.collapsePointsOnCellEdges(res.coordinates, patch, singularIds, + opts_.contourAlignmentAngle, edgeIds); + sT_.collapsePointsOnFeatureEdges(res.coordinates, patch, singularIds, + incidentElements, edgeIds); + } + meshTools::checkNoCellsAreCrossed(res); - std::for_each( -#ifdef TESSELLATOR_EXECUTION_POLICIES - std::execution::par, -#endif - patchs.begin(), patchs.end(), [&](auto& p) { - sT_.collapseInteriorPointsToBound(res.coordinates, p); - }); + for (const auto& patch : patchs) { + sT_.collapsePointsOnCellFaces(res.coordinates, patch, singularIds, faceIds); + sT_.collapsePointsOnFeatureEdges(res.coordinates, patch, singularIds, + incidentElements, faceIds); + } + meshTools::checkNoCellsAreCrossed(res); + + for (const auto& patch : patchs) { + sT_.collapsePointsOnFeatureEdges(res.coordinates, patch, singularIds, + incidentElements, interiorIds); + sT_.collapseInteriorPointsToBound(res.coordinates, patch, interiorIds); + } + meshTools::checkNoCellsAreCrossed(res); } @@ -95,4 +168,4 @@ Smoother::Smoother(const Mesh& mesh, const SmootherOptions& opts) : } } -} \ No newline at end of file +} diff --git a/src/core/SmootherTools.cpp b/src/core/SmootherTools.cpp index fd38246..dc46199 100644 --- a/src/core/SmootherTools.cpp +++ b/src/core/SmootherTools.cpp @@ -57,10 +57,66 @@ void SmootherTools::updateCoordinates( } } +bool SmootherTools::moveWouldCrossGrid( + const CoordinateId& id, + const Coordinate& destination, + const Coordinates& coordinates, + const IncidentElements& incidentElements) const +{ + const auto incident = incidentElements.find(id); + if (incident == incidentElements.end()) { + return false; + } + + for (const auto* element : incident->second) { + std::set commonCells; + bool firstVertex = true; + for (const auto vertexId : element->vertices) { + const auto& coordinate = vertexId == id ? destination : coordinates[vertexId]; + const auto touchingCells = getTouchingCells(coordinate); + if (firstVertex) { + commonCells = touchingCells; + firstVertex = false; + } + else { + for (auto cell = commonCells.begin(); cell != commonCells.end();) { + if (touchingCells.count(*cell) == 0) { + cell = commonCells.erase(cell); + } + else { + ++cell; + } + } + } + if (commonCells.empty()) { + return true; + } + } + } + return false; +} + +void SmootherTools::collapsePointsOnFeatureEdges( + Coordinates& coords, + const ElementsView& patch, + const SingularIds& singularIds, + const IdSet& movableIds) +{ + IncidentElements incidentElements; + for (const auto* element : patch) { + for (const auto id : element->vertices) { + incidentElements[id].push_back(element); + } + } + collapsePointsOnFeatureEdges(coords, patch, singularIds, incidentElements, movableIds); +} + void SmootherTools::collapsePointsOnFeatureEdges( Coordinates& coords, const ElementsView& patch, - const SingularIds& singularIds) + const SingularIds& singularIds, + const IncidentElements& incidentElements, + const IdSet& movableIds) { CoordGraph Point = CoordGraph(patch); CoordGraph edges = Point.getBoundaryGraph().intersect(singularIds.featureIds()); @@ -98,11 +154,31 @@ void SmootherTools::collapsePointsOnFeatureEdges( std::map toMove; for (auto const& i : validInterior) { + if (!movableIds.empty() && movableIds.count(i) == 0) { + continue; + } if (isRelativeInCellCorner(coords[i])) { continue; } - Coordinate closest = closestByDistance(coords, i, Point.getClosestVerticesInSet(i, validExterior)); + auto candidates = Point.getClosestVerticesInSet(i, validExterior); + const auto sourceCells = getTouchingCells(coords[i]); + for (auto candidate = candidates.begin(); candidate != candidates.end();) { + const auto candidateCells = getTouchingCells(coords[*candidate]); + const bool sharesCell = std::any_of( + candidateCells.begin(), candidateCells.end(), + [&](const Cell& cell) { return sourceCells.count(cell) != 0; }); + if (sharesCell) { + ++candidate; + } + else { + candidate = candidates.erase(candidate); + } + } + if (candidates.empty()) { + continue; + } + Coordinate closest = closestByDistance(coords, i, candidates); if (isRelativeInCellFace(coords[i]) && !areCoordOnSameFace(coords[i], closest)) { continue; } @@ -112,7 +188,12 @@ void SmootherTools::collapsePointsOnFeatureEdges( toMove[i] = closest; } - updateCoordinates(coords, toMove); + std::lock_guard lock(writingCoordinates_); + for (const auto& move : toMove) { + if (!moveWouldCrossGrid(move.first, move.second, coords, incidentElements)) { + coords[move.first] = move.second; + } + } } Coordinate SmootherTools::closestByDistance( @@ -169,7 +250,8 @@ void SmootherTools::collapsePointsOnCellEdges( Coordinates& coords, const ElementsView& patch, const SingularIds& singularIds, - double alignmentAngle) + double alignmentAngle, + const IdSet& movableIds) { { IdSet vertices = CoordGraph(patch).getVertices(); @@ -199,7 +281,10 @@ void SmootherTools::collapsePointsOnCellEdges( IdSet interiorValid = intersectWithIdSet(interior, protectedIds); - IdSet movable = classifyIds(interior, [&](auto i) {return !protectedIds.count(i); }).first; + IdSet movable = classifyIds(interior, [&](auto i) { + return !protectedIds.count(i) + && (movableIds.empty() || movableIds.count(i) != 0); + }).first; IdSet validIds = mergeIds(cG.getExterior(), interiorValid); std::map toMove; @@ -207,7 +292,11 @@ void SmootherTools::collapsePointsOnCellEdges( if (isRelativeInCellCorner(coords[i])) { continue; } - Coordinate closest = closestByDistance(coords, i, cG.getClosestVerticesInSet(i, validIds)); + const auto candidates = cG.getClosestVerticesInSet(i, validIds); + if (candidates.empty()) { + continue; + } + Coordinate closest = closestByDistance(coords, i, candidates); if (isRelativeInCellFace(coords[i]) && !areCoordOnSameFace(coords[i], closest)) { continue; } @@ -241,7 +330,8 @@ IdSet SmootherTools::getClosestValidByDistanceInCycle( void SmootherTools::collapsePointsOnCellFaces( Coordinates& coords, const ElementsView& patch, - const SingularIds& sIds) + const SingularIds& sIds, + const IdSet& movableIds) { std::map toMove; std::map > cyclesToValidOrOnFace; @@ -260,6 +350,9 @@ void SmootherTools::collapsePointsOnCellFaces( const IdSet& valid = kv.second.first; const IdSet& onCellFace = kv.second.second; for (auto const& id : onCellFace) { + if (!movableIds.empty() && movableIds.count(id) == 0) { + continue; + } try { IdSet candidates = getClosestValidByDistanceInCycle(id, cycle, valid); if (!candidates.empty()) { @@ -426,7 +519,7 @@ CoordGraph::Path SmootherTools::pathFromIdToAnyTarget( } const std::size_t i = it - cycle.begin(); - CoordGraph::Path res(startId); + CoordGraph::Path res{startId}; for (std::size_t d = 1; d < cycle.size(); d++) { CoordinateId idTest; if (forward) { @@ -485,13 +578,17 @@ Coordinates SmootherTools::collapsePointsOnContour( void SmootherTools::collapseInteriorPointsToBound( Coordinates& coords, - const ElementsView& patch) + const ElementsView& patch, + const IdSet& movableIds) { IdSet bound, interior; std::tie(bound, interior) = CoordGraph(patch).getBoundAndInteriorVertices(); std::map toMove; for (auto const& vI : interior) { + if (!movableIds.empty() && movableIds.count(vI) == 0) { + continue; + } toMove[vI] = closestByDistance(coords, vI, bound); } @@ -573,4 +670,4 @@ void SmootherTools::reorientSingleElement( } } -} \ No newline at end of file +} diff --git a/src/core/SmootherTools.h b/src/core/SmootherTools.h index 973de4f..3ab957b 100644 --- a/src/core/SmootherTools.h +++ b/src/core/SmootherTools.h @@ -13,6 +13,8 @@ namespace core { class SmootherTools : public utils::GridTools { public: + using IncidentElements = std::map; + class SingularIds { public: SingularIds(const IdSet& featureIds, const IdSet& contourIds, const IdSet& cornerIds) : @@ -39,7 +41,15 @@ class SmootherTools : public utils::GridTools { void collapsePointsOnFeatureEdges( Coordinates& res, const ElementsView& patch, - const SingularIds& singularIds); + const SingularIds& singularIds, + const IdSet& movableIds = {}); + + void collapsePointsOnFeatureEdges( + Coordinates& res, + const ElementsView& patch, + const SingularIds& singularIds, + const IncidentElements& incidentElements, + const IdSet& movableIds = {}); Coordinates collapsePointsOnContour( const Elements& elems, @@ -50,12 +60,14 @@ class SmootherTools : public utils::GridTools { Coordinates& res, const ElementsView& patch, const SingularIds& singularIds , - double alignmentAngle); + double alignmentAngle, + const IdSet& movableIds = {}); void collapsePointsOnCellFaces( Coordinates& res, const ElementsView& patch, - const SingularIds&); + const SingularIds&, + const IdSet& movableIds = {}); void remeshBoundary( Elements& es, @@ -70,7 +82,8 @@ class SmootherTools : public utils::GridTools { void collapseInteriorPointsToBound( Coordinates& coords, - const ElementsView& patch); + const ElementsView& patch, + const IdSet& movableIds = {}); void remeshElementsToOneInteriorPoint( Elements& es, @@ -83,11 +96,19 @@ class SmootherTools : public utils::GridTools { const ElementsView& patch); private: + friend class SmootherToolsTestAccess; + std::mutex writingCoordinates_; std::mutex writingElements_; void updateCoordinates(Coordinates& res, std::map toMove); + bool moveWouldCrossGrid( + const CoordinateId& id, + const Coordinate& destination, + const Coordinates& coordinates, + const IncidentElements& incidentElements) const; + static CoordinateId getClosestEndOfPaths( const std::vector& paths); @@ -136,4 +157,4 @@ class SmootherTools : public utils::GridTools { }; } -} \ No newline at end of file +} diff --git a/src/core/Staircaser.cpp b/src/core/Staircaser.cpp index 40c9d00..f51f2af 100644 --- a/src/core/Staircaser.cpp +++ b/src/core/Staircaser.cpp @@ -448,7 +448,7 @@ bool Staircaser::isAllCoordinatesOnTheSameCellBoundary(const Element& newElement return !commonStructuredCells.empty(); } -std::set Staircaser::getElementsConvertedInLines(Cell cell, std::map> cellElemMap) { +std::set Staircaser::getElementsConvertedInLines(const Cell cell, const std::map>& cellElemMap) { std::set elementsConvertedInLines; for (const auto& e : cellElemMap.at(cell)) { @@ -492,7 +492,7 @@ std::set Staircaser::getElementsConvertedInLines(Cell cell, std::map elementsConvertedInLines, Group& meshGroup, const Cell cell, std::map> cellElemMap, std::vector> toRemove) { +void Staircaser::splitLinesWithNeighborTriangle(const size_t groupIndex, std::set elementsConvertedInLines, Group& meshGroup, const Cell cell, const std::map>& cellElemMap, std::vector> toRemove) { for (const auto& e: elementsConvertedInLines) { bool processed = false; for (const auto& otherElementsInCell: cellElemMap.at(cell)) { diff --git a/src/core/Staircaser.h b/src/core/Staircaser.h index a0a94e4..c836822 100644 --- a/src/core/Staircaser.h +++ b/src/core/Staircaser.h @@ -78,10 +78,12 @@ class Staircaser : public utils::GridTools { void fillGaps(const RelativePairSet boundaryCoordinatePairs); std::pair obtainNewIndexForElement(const Element& e, const std::set& cellsToStructure, CoordinateMap& coordinateMap); bool isAllCoordinatesOnTheSameCellBoundary(const Element& newElement, const std::set& cellsToStructure); - std::set getElementsConvertedInLines(const Cell cell, std::map> cellElemMap); - void splitLinesWithNeighborTriangle(const size_t groupIndex, std::set elementsConvertedInLines, Group& meshGroup, const Cell cell, std::map> cellElemMap, std::vector> toRemove); + std::set getElementsConvertedInLines(const Cell cell, const std::map>& cellElemMap); + void splitLinesWithNeighborTriangle(const size_t groupIndex, std::set elementsConvertedInLines, Group& meshGroup, const Cell cell, const std::map>& cellElemMap, std::vector> toRemove); }; + + } } \ No newline at end of file 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 dfdd88c..b48f2ba 100644 --- a/src/meshers/CMakeLists.txt +++ b/src/meshers/CMakeLists.txt @@ -6,10 +6,12 @@ add_library(tessellator-meshers "OffgridMesher.cpp" "ConformalMesher.cpp" ) -target_link_libraries(tessellator-meshers tessellator-core tessellator-utils) +target_link_libraries(tessellator-meshers + tessellator-core + 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/ConformalMesher.h b/src/meshers/ConformalMesher.h index 8fe5540..26a9858 100644 --- a/src/meshers/ConformalMesher.h +++ b/src/meshers/ConformalMesher.h @@ -17,6 +17,8 @@ class ConformalMesher : public MesherBase { virtual ~ConformalMesher() = default; Mesh mesh() const; + + const ConformalMesherOptions & getOptions() const { return opts_; } static std::set findNonConformalCells(const Mesh& mesh); static std::set cellsWithMoreThanAVertexInsideEdge(const Mesh& mesh); diff --git a/src/meshers/ConformalMesherOptions.h b/src/meshers/ConformalMesherOptions.h index 86f31d2..5a38c18 100644 --- a/src/meshers/ConformalMesherOptions.h +++ b/src/meshers/ConformalMesherOptions.h @@ -1,14 +1,15 @@ #pragma once #include "types/Mesh.h" +#include "MesherBaseOptions.h" #include "core/SnapperOptions.h" namespace meshlib::meshers { -class ConformalMesherOptions { +class ConformalMesherOptions : public MesherBaseOptions { public: core::SnapperOptions snapperOptions; - std::set volumeGroups{}; + // std::set volumeGroups{}; }; } 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 a688c29..d05edd2 100644 --- a/src/meshers/MesherBase.h +++ b/src/meshers/MesherBase.h @@ -1,6 +1,7 @@ #pragma once #include "types/Mesh.h" +#include "MesherBaseOptions.h" namespace meshlib { namespace meshers { @@ -10,6 +11,7 @@ class MesherBase { MesherBase(const Mesh& in); virtual ~MesherBase() = default; virtual Mesh mesh() const = 0; + const MesherBaseOptions & getOptions() const { return opts_; } protected: virtual void process(Mesh&) const = 0; @@ -19,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); @@ -29,7 +32,10 @@ class MesherBase { Grid originalGrid_; Grid enlargedGrid_; + MesherBaseOptions opts_; + + }; } -} \ No newline at end of file +} diff --git a/src/meshers/MesherBaseOptions.h b/src/meshers/MesherBaseOptions.h new file mode 100644 index 0000000..55138c1 --- /dev/null +++ b/src/meshers/MesherBaseOptions.h @@ -0,0 +1,13 @@ +#pragma once + +#include "types/Mesh.h" +#include "core/SnapperOptions.h" + +namespace meshlib::meshers { + +class MesherBaseOptions { +public: + std::set volumeGroups{}; +}; + +} diff --git a/src/meshers/StaircaseMesher.cpp b/src/meshers/StaircaseMesher.cpp index 7ad6458..fadfc34 100644 --- a/src/meshers/StaircaseMesher.cpp +++ b/src/meshers/StaircaseMesher.cpp @@ -6,6 +6,9 @@ #include "core/Slicer.h" #include "core/Collapser.h" #include "core/Staircaser.h" +#include "core/Compressor.h" +#include "core/VolumeFiller.h" +#include "core/VolumeShellExtractor.h" #include "utils/RedundancyCleaner.h" #include "utils/MeshTools.h" @@ -17,17 +20,35 @@ using namespace utils; using namespace core; using namespace meshTools; -StaircaseMesher::StaircaseMesher(const Mesh& inputMesh, int decimalPlacesInCollapser) : +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) + 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) @@ -37,13 +58,40 @@ 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; } @@ -71,6 +119,24 @@ void StaircaseMesher::process(Mesh& mesh) const logNumberOfQuads(countMeshElementsIf(mesh, isQuad)); logNumberOfLines(countMeshElementsIf(mesh, isLine)); + + if (compress) { + log("Compressing surfaces.", 1); + std::size_t beforeQuads = countMeshElementsIf(mesh, isQuad); + std::size_t merged = Compressor::compressSurfacesInMesh(mesh); + std::size_t afterQuads = countMeshElementsIf(mesh, isQuad); + log("Compressed " + std::to_string(beforeQuads) + + " -> " + std::to_string(afterQuads) + + " quads (merged " + std::to_string(merged) + " surfaces)", 1); + + log("Compressing lines.", 1); + std::size_t beforeLines = countMeshElementsIf(mesh, isLine); + merged = Compressor::compressLinesInMesh(mesh, dimensions); + std::size_t afterLines = countMeshElementsIf(mesh, isLine); + log("Compressed " + std::to_string(beforeLines) + + " -> " + std::to_string(afterLines) + + " lines (merged " + std::to_string(merged) + " segments)", 1); + } log("Recovering original grid size.", 1); reduceGrid(mesh, originalGrid_); @@ -81,6 +147,8 @@ void StaircaseMesher::process(Mesh& mesh) const logNumberOfQuads(countMeshElementsIf(mesh, isQuad)); logNumberOfLines(countMeshElementsIf(mesh, isLine)); + copyGroupNames(mesh, groupNames); + } @@ -89,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 bbaffa1..d768399 100644 --- a/src/meshers/StaircaseMesher.h +++ b/src/meshers/StaircaseMesher.h @@ -2,22 +2,27 @@ #include "types/Mesh.h" #include "MesherBase.h" +#include "StaircaseMesherOptions.h" namespace meshlib::meshers { class StaircaseMesher : public MesherBase { public: - StaircaseMesher(const Mesh& in, int decimalPlacesInCollapser = 4); + StaircaseMesher(const Mesh& in, int decimalPlacesInCollapser = 4, StaircaseMesherOptions opts = StaircaseMesherOptions()); virtual ~StaircaseMesher() = default; Mesh mesh() const; + const StaircaseMesherOptions & getOptions() const { return opts_; } private: 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 new file mode 100644 index 0000000..4bac022 --- /dev/null +++ b/src/meshers/StaircaseMesherOptions.h @@ -0,0 +1,13 @@ +#pragma once + +#include "MesherBaseOptions.h" + +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 bb0d41a..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; @@ -134,10 +139,15 @@ typedef std::size_t ElementId; typedef std::vector Elements; struct Group { + std::string name; std::vector elements; + Group() = default; + Group(const std::vector& elems) : elements(elems) {} + Group(const std::string& n, const std::vector& elems) : name(n), elements(elems) {} + bool operator==(const Group& rhs) const { - return elements == rhs.elements; + return name == rhs.name && elements == rhs.elements; } std::map> buildCoordToElemMap() const { @@ -155,6 +165,7 @@ struct Group { friend class boost::serialization::access; template void serialize(Archive& ar, const unsigned int version) { + ar& name; ar& elements; } }; @@ -214,4 +225,3 @@ struct Mesh { }; } - diff --git a/src/utils/CoordGraph.cpp b/src/utils/CoordGraph.cpp index e4242f6..0e5e0c3 100644 --- a/src/utils/CoordGraph.cpp +++ b/src/utils/CoordGraph.cpp @@ -81,6 +81,10 @@ CoordGraph::CoordGraph(const Elements& elems) CoordGraph::CoordGraph(const ElementsView& es) { for (auto const& e : es) { + if (e->vertices.size() == 1) { + addVertex(e->vertices.front()); + continue; + } for (std::size_t i = 0; i < e->vertices.size(); i++) { this->addEdge( e->vertices[i], @@ -95,6 +99,10 @@ CoordGraph::CoordGraph(const ElementsView& es) CoordGraph::CoordGraph(const Paths& paths) { for (const auto& p : paths) { + if (p.size() == 1) { + addVertex(p.front()); + continue; + } for (std::size_t i = 0; i < p.size(); i++) { this->addEdge( p[i], @@ -221,7 +229,7 @@ IdSet CoordGraph::getClosestVerticesInSet( } auto path = findShortestPath(vI, vB); if (path.empty()) { - throw std::runtime_error("Can not find path to point in set."); + continue; } paths.insert(path); } @@ -586,4 +594,4 @@ std::vector CoordGraph::findCycles() const } } -} \ No newline at end of file +} diff --git a/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 b2a25f4..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()); @@ -418,4 +440,54 @@ bool isAClosedTopology(const Elements& es) return CoordGraph(es).getBoundaryGraph().getVertices().size() == 0; } +Mesh extractGroupsByName(const Mesh& mesh, const std::vector& groupNames) +{ + Mesh result; + result.grid = mesh.grid; + + std::map coordRemap; + std::map> groupCoordIds; + + for (const auto& groupName : groupNames) { + auto it = std::find_if(mesh.groups.begin(), mesh.groups.end(), + [&groupName](const Group& g) { return g.name == groupName; }); + + if (it != mesh.groups.end()) { + result.groups.push_back(*it); + + for (const auto& elem : it->elements) { + for (const auto& vId : elem.vertices) { + std::string coordKey = groupName + "_" + std::to_string(vId); + if (coordRemap.find(coordKey) == coordRemap.end()) { + CoordinateId newVId = result.coordinates.size(); + result.coordinates.push_back(mesh.coordinates[vId]); + coordRemap[coordKey] = newVId; + groupCoordIds[groupName].push_back(vId); + } + } + } + } else { + result.groups.push_back(Group{groupName, {}}); + } + } + + for (auto& group : result.groups) { + std::string groupName = group.name; + std::map elemRemap; + + for (const auto& oldVId : groupCoordIds[groupName]) { + std::string coordKey = groupName + "_" + std::to_string(oldVId); + elemRemap[oldVId] = coordRemap[coordKey]; + } + + for (auto& elem : group.elements) { + for (auto& vId : elem.vertices) { + vId = elemRemap[vId]; + } + } + } + + return result; +} + } \ No newline at end of file diff --git a/src/utils/MeshTools.h b/src/utils/MeshTools.h index 688c0cb..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); @@ -50,4 +53,6 @@ void mergeMeshAsNewGroup(Mesh& lMesh, const Mesh& iMesh); bool isAClosedTopology(const Elements& es); -} \ No newline at end of file +Mesh extractGroupsByName(const Mesh& mesh, const std::vector& groupNames); + +} diff --git a/src/utils/Types.h b/src/utils/Types.h index 471ca81..e5f1668 100644 --- a/src/utils/Types.h +++ b/src/utils/Types.h @@ -61,5 +61,17 @@ using HexIds = std::array; using UpdateMap = std::array, 2>, 2>; +// Compressor types +using Sign = int; +using SignedAxis = std::pair; +using GridLine = std::array; + +using PlanePoint = std::array; +using PlaneLinel = std::pair; +using PlaneSurfel = PlanePoint; +using PlaneSurface = std::pair; +using Contour = std::vector; +using CrossLine = std::pair>; + } diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index cea5caa..3cc1202 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -3,7 +3,6 @@ message(STATUS "Creating build system for tessellator-tests") find_package(GTest CONFIG REQUIRED) include_directories(${GTEST_INCLUDE_DIRS}) - include_directories( ${PROJECT_SOURCE_DIR}/src/ ${PROJECT_SOURCE_DIR}/src/app/ @@ -15,14 +14,15 @@ include_directories( ) add_executable(tessellator_tests - "app/launcherTest.cpp" - "app/vtkIOTest.cpp" "core/CollapserTest.cpp" + "core/CompressorTest.cpp" "core/SlicerTest.cpp" - "core/SnapperTest.cpp" "core/SmootherTest.cpp" "core/SmootherToolsTest.cpp" - "core/StaircaserTest.cpp" + "core/SnapperTest.cpp" + "core/StaircaserTest.cpp" + "core/VolumeFillerTest.cpp" + "core/VolumeShellExtractorTest.cpp" "types/MeshTest.cpp" "utils/ConvexHullTest.cpp" "utils/CoordGraphTest.cpp" @@ -31,18 +31,25 @@ add_executable(tessellator_tests "utils/GridToolsTest.cpp" "utils/MeshToolsTest.cpp" "utils/RedundancyCleanerTest.cpp" - "meshers/StaircaseMesherTest.cpp" - "meshers/OffgridMesherTest.cpp" "meshers/ConformalMesherTest.cpp" + "meshers/OffgridMesherTest.cpp" + "meshers/StaircaseMesherTest.cpp" ) target_link_libraries(tessellator_tests tessellator-meshers - tessellator-app GTest::gtest GTest::gtest_main ) +if(TESSELLATOR_LOAD_APP) + target_sources(tessellator_tests PRIVATE + "app/launcherTest.cpp" + "app/vtkIOTest.cpp" + ) + target_link_libraries(tessellator_tests tessellator-app) +endif() + if (TESSELLATOR_ENABLE_CGAL) include_directories( ${PROJECT_SOURCE_DIR}/src/cgal/ @@ -62,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 bda4612..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; @@ -1269,4 +1270,49 @@ static Mesh buildProblematicTriMesh2() } } + +static CoordinateId findOrAddCoord(Mesh& mesh, const std::array& gridIdx) { + double pos[3]; + pos[0] = mesh.grid[0][gridIdx[0]]; + pos[1] = mesh.grid[1][gridIdx[1]]; + pos[2] = mesh.grid[2][gridIdx[2]]; + + for (CoordinateId i = 0; i < static_cast(mesh.coordinates.size()); ++i) { + if (mesh.coordinates[i](0) == pos[0] && + mesh.coordinates[i](1) == pos[1] && + mesh.coordinates[i](2) == pos[2]) { + return i; + } + } + mesh.coordinates.push_back(Coordinate({pos[0], pos[1], pos[2]})); + return static_cast(mesh.coordinates.size() - 1); +} + +// Helper to add a quad (as a surface with four vertices) to a mesh +static void addQuad(Mesh& mesh, const std::array& v0, const std::array& v1, + const std::array& v2, const std::array& v3, GroupId groupId = 0) { + if (mesh.groups.size() <= groupId) { + mesh.groups.resize(groupId + 1); + } + + CoordinateId c0 = findOrAddCoord(mesh, v0); + CoordinateId c1 = findOrAddCoord(mesh, v1); + CoordinateId c2 = findOrAddCoord(mesh, v2); + CoordinateId c3 = findOrAddCoord(mesh, v3); + + mesh.groups[groupId].elements.push_back(Element({c0, c1, c2, c3}, Element::Type::Surface)); +} + +// Helper to add a line (as a line with two vertices) to a mesh +static void addLine(Mesh& mesh, const std::array& v0, const std::array& v1, GroupId groupId = 0) { + if (mesh.groups.size() <= groupId) { + mesh.groups.resize(groupId + 1); + } + + CoordinateId c0 = findOrAddCoord(mesh, v0); + CoordinateId c1 = findOrAddCoord(mesh, v1); + + mesh.groups[groupId].elements.push_back(Element({c0, c1}, Element::Type::Line)); +} + } \ No newline at end of file diff --git a/test/app/launcherTest.cpp b/test/app/launcherTest.cpp index 3942428..2cf01ab 100644 --- a/test/app/launcherTest.cpp +++ b/test/app/launcherTest.cpp @@ -2,6 +2,8 @@ #include "app/launcher.h" #include "types/Mesh.h" +#include "meshers/StaircaseMesher.h" + #include #include @@ -15,7 +17,7 @@ TEST_F(LauncherTest, prints_help) { int ac = 2; const char* av[] = { NULL, "-h" }; - EXPECT_EQ(meshlib::app::launcher(ac, av), EXIT_SUCCESS); + EXPECT_EQ(launcher(ac, av), EXIT_SUCCESS); } TEST_F(LauncherTest, parse_rectilinear_grid) @@ -28,7 +30,7 @@ TEST_F(LauncherTest, parse_rectilinear_grid) i >> j; } - meshlib::Grid grid = meshlib::app::parseGridFromJSON(j["grid"]); + meshlib::Grid grid = parseGridFromJSON(j["grid"]); meshlib::Grid expectedGrid({ std::vector{600, 603.25}, @@ -49,12 +51,108 @@ TEST_F(LauncherTest, parse_rectilinear_grid) } } +TEST_F(LauncherTest, builds_staircased_mesher_default) +{ + meshlib::Mesh meshMock; + meshMock.grid = { + std::vector{0, 1}, + std::vector{0, 1}, + std::vector{0, 1} + }; + std::string fileName = "testData/cases/longPolyline/longPolyline_legacy.tessellator.json"; + nlohmann::json j; + { + std::ifstream i(fileName); + i >> j; + } + + ObjectDefinition objDef; + auto mesher = buildMesher(meshMock, j, objDef); + + EXPECT_NO_THROW(auto staircaseMesher = dynamic_cast(*mesher)); + const auto & options = dynamic_cast(*mesher).getOptions(); + 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) +{ + meshlib::Mesh meshMock; + meshMock.grid = { + std::vector{0, 1}, + std::vector{0, 1}, + std::vector{0, 1} + }; + std::string fileName = "testData/cases/longPolyline/longPolyline.tessellator.json"; + nlohmann::json j; + { + std::ifstream i(fileName); + i >> j; + } + + ObjectDefinition objDef; + auto mesher = buildMesher(meshMock, j, objDef); + + EXPECT_NO_THROW(auto staircaseMesher = dynamic_cast(*mesher)); + const auto & options = dynamic_cast(*mesher).getOptions(); + EXPECT_EQ(options.volumeGroups.size(), 0); + EXPECT_EQ(options.compress, false); +} + +TEST_F(LauncherTest, builds_staircased_mesher_with_compression) +{ + meshlib::Mesh meshMock; + meshMock.grid = { + std::vector{0, 1}, + std::vector{0, 1}, + std::vector{0, 1} + }; + std::string fileName = "testData/cases/longPolyline/longPolyline_compression.tessellator.json"; + nlohmann::json j; + { + std::ifstream i(fileName); + i >> j; + } + + ObjectDefinition objDef; + auto mesher = buildMesher(meshMock, j, objDef); + + EXPECT_NO_THROW(auto staircaseMesher = dynamic_cast(*mesher)); + const auto & options = dynamic_cast(*mesher).getOptions(); + EXPECT_EQ(options.volumeGroups.size(), 0); + EXPECT_EQ(options.compress, true); +} + TEST_F(LauncherTest, launches_alhambra_case) { int ac = 3; const char* av[] = { NULL, "-i", "testData/cases/alhambra/alhambra.tessellator.json"}; int exitCode; - EXPECT_NO_THROW(exitCode = meshlib::app::launcher(ac, av)); + EXPECT_NO_THROW(exitCode = launcher(ac, av)); EXPECT_EQ(exitCode, EXIT_SUCCESS); } @@ -63,7 +161,7 @@ TEST_F(LauncherTest, launches_conformal_alhambra_case) int ac = 3; const char* av[] = { NULL, "-i", "testData/cases/alhambra/alhambra.conformal.tessellator.json"}; int exitCode; - EXPECT_NO_THROW(exitCode = meshlib::app::launcher(ac, av)); + EXPECT_NO_THROW(exitCode = launcher(ac, av)); EXPECT_EQ(exitCode, EXIT_SUCCESS); } @@ -73,7 +171,16 @@ TEST_F(LauncherTest, launches_sphere_case) int ac = 3; const char* av[] = { NULL, "-i", "testData/cases/sphere/sphere.tessellator.json"}; int exitCode; - EXPECT_NO_THROW(exitCode = meshlib::app::launcher(ac, av)); + EXPECT_NO_THROW(exitCode = launcher(ac, av)); + EXPECT_EQ(exitCode, EXIT_SUCCESS); +} + +TEST_F(LauncherTest, launches_closed_sphere_case) +{ + int ac = 3; + const char* av[] = { NULL, "-i", "testData/cases/sphere/closed_sphere.tessellator.json"}; + int exitCode; + EXPECT_NO_THROW(exitCode = launcher(ac, av)); EXPECT_EQ(exitCode, EXIT_SUCCESS); } @@ -82,7 +189,7 @@ TEST_F(LauncherTest, launches_conformal_sphere_case) int ac = 3; const char* av[] = { NULL, "-i", "testData/cases/sphere/sphere.conformal.tessellator.json"}; int exitCode; - EXPECT_NO_THROW(exitCode = meshlib::app::launcher(ac, av)); + EXPECT_NO_THROW(exitCode = launcher(ac, av)); EXPECT_EQ(exitCode, EXIT_SUCCESS); } @@ -91,7 +198,7 @@ TEST_F(LauncherTest, launches_conformal_thinCylinder_case) int ac = 3; const char* av[] = { NULL, "-i", "testData/cases/thinCylinder/thinCylinder.conformal.tessellator.json"}; int exitCode; - EXPECT_NO_THROW(exitCode = meshlib::app::launcher(ac, av)); + EXPECT_NO_THROW(exitCode = launcher(ac, av)); EXPECT_EQ(exitCode, EXIT_SUCCESS); } @@ -100,7 +207,7 @@ TEST_F(LauncherTest, launches_thinCylinder_case) int ac = 3; const char* av[] = { NULL, "-i", "testData/cases/thinCylinder/thinCylinder.tessellator.json"}; int exitCode; - EXPECT_NO_THROW(exitCode = meshlib::app::launcher(ac, av)); + EXPECT_NO_THROW(exitCode = launcher(ac, av)); EXPECT_EQ(exitCode, EXIT_SUCCESS); } @@ -109,16 +216,16 @@ TEST_F(LauncherTest, launches_cone_case) int ac = 3; const char* av[] = { NULL, "-i", "testData/cases/cone/cone.tessellator.json" }; int exitCode; - EXPECT_NO_THROW(exitCode = meshlib::app::launcher(ac, av)); + EXPECT_NO_THROW(exitCode = launcher(ac, av)); EXPECT_EQ(exitCode, EXIT_SUCCESS); } TEST_F(LauncherTest, launches_long_polyline_case) { int ac = 3; - const char* av[] = { NULL, "-i", "testData/cases/longPolyline/longPolyline.tessellator.json" }; + const char* av[] = { NULL, "-i", "testData/cases/longPolyline/longPolyline_legacy.tessellator.json" }; int exitCode; - EXPECT_NO_THROW(exitCode = meshlib::app::launcher(ac, av)); + EXPECT_NO_THROW(exitCode = launcher(ac, av)); EXPECT_EQ(exitCode, EXIT_SUCCESS); } @@ -127,7 +234,147 @@ TEST_F(LauncherTest, launches_conformal_cone_case) int ac = 3; const char* av[] = { NULL, "-i", "testData/cases/cone/cone.conformal.tessellator.json"}; int exitCode; - EXPECT_NO_THROW(exitCode = meshlib::app::launcher(ac, av)); + EXPECT_NO_THROW(exitCode = launcher(ac, av)); EXPECT_EQ(exitCode, EXIT_SUCCESS); } +TEST_F(LauncherTest, readObjectsFromJSON_basic) +{ + std::string fileName = "testData/cases/multiObject/basic.tessellator.json"; + nlohmann::json j; + { + std::ifstream i(fileName); + i >> j; + } + auto objects = readObjectsFromJSON(j); + EXPECT_EQ(objects.size(), 2); + EXPECT_EQ(objects[0].filename, "sphere.stl"); + EXPECT_EQ(objects[0].group, "sphere_group"); + EXPECT_TRUE(objects[0].isVolume); + 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].mesherOverride.has_value()); +} + +TEST_F(LauncherTest, readObjectsFromJSON_mixedMesher) +{ + std::string fileName = "testData/cases/multiObject/mixedMesher.tessellator.json"; + nlohmann::json j; + { + std::ifstream i(fileName); + i >> j; + } + auto objects = readObjectsFromJSON(j); + EXPECT_EQ(objects.size(), 3); + + EXPECT_EQ(objects[0].filename, "sphere.stl"); + EXPECT_FALSE(objects[0].isVolume); + EXPECT_TRUE(objects[0].mesherOverride.has_value()); + EXPECT_EQ(objects[0].mesherOverride.value()["type"], "staircase"); + EXPECT_TRUE(objects[0].mesherOverride.value().contains("options")); + EXPECT_TRUE(objects[0].mesherOverride.value()["options"].contains("compress")); + EXPECT_TRUE(objects[0].mesherOverride.value()["options"]["compress"]); + + EXPECT_EQ(objects[1].filename, "cone.stl"); + EXPECT_FALSE(objects[1].isVolume); + EXPECT_TRUE(objects[1].mesherOverride.has_value()); + EXPECT_EQ(objects[1].mesherOverride.value()["type"], "conformal"); + + EXPECT_EQ(objects[2].filename, "cone.stl"); + EXPECT_FALSE(objects[2].isVolume); + EXPECT_FALSE(objects[2].mesherOverride.has_value()); +} + +TEST_F(LauncherTest, readObjectsFromJSON_singleObject) +{ + std::string fileName = "testData/cases/multiObject/singleObject.tessellator.json"; + nlohmann::json j; + { + std::ifstream i(fileName); + i >> j; + } + auto objects = readObjectsFromJSON(j); + EXPECT_EQ(objects.size(), 1); + EXPECT_EQ(objects[0].filename, "sphere.stl"); + EXPECT_EQ(objects[0].group, "default_group"); + EXPECT_FALSE(objects[0].isVolume); +} + +TEST_F(LauncherTest, readObjectsFromJSON_legacyFormat) +{ + std::string fileName = "testData/cases/sphere/closed_sphere.tessellator.json"; + nlohmann::json j; + { + std::ifstream i(fileName); + i >> j; + } + auto objects = readObjectsFromJSON(j); + EXPECT_EQ(objects.size(), 1); + EXPECT_EQ(objects[0].filename, "sphere.stl"); + EXPECT_EQ(objects[0].group, "sphere"); + EXPECT_TRUE(objects[0].isVolume); +} + +TEST_F(LauncherTest, builds_staircased_mesher_with_override) +{ + meshlib::Mesh meshMock; + meshMock.grid = { + std::vector{0, 1}, + std::vector{0, 1}, + std::vector{0, 1} + }; + + std::string fileName = "testData/cases/multiObject/mixedMesher.tessellator.json"; + nlohmann::json j; + { + std::ifstream i(fileName); + i >> j; + } + + auto objects = readObjectsFromJSON(j); + + auto mesher = meshlib::app::buildMesher(meshMock, j, objects[0]); + + EXPECT_NO_THROW(auto staircaseMesher = dynamic_cast(*mesher)); + const auto & options = dynamic_cast(*mesher).getOptions(); + EXPECT_EQ(options.volumeGroups.size(), 0); + EXPECT_EQ(options.compress, true); +} + +TEST_F(LauncherTest, launches_multiObject_basic) +{ + int ac = 3; + const char* av[] = { NULL, "-i", "testData/cases/multiObject/basic.tessellator.json"}; + int exitCode; + EXPECT_NO_THROW(exitCode = launcher(ac, av)); + EXPECT_EQ(exitCode, EXIT_SUCCESS); +} + +TEST_F(LauncherTest, launches_multiObject_mixedMesher) +{ + int ac = 3; + const char* av[] = { NULL, "-i", "testData/cases/multiObject/mixedMesher.tessellator.json"}; + int exitCode; + EXPECT_NO_THROW(exitCode = launcher(ac, av)); + EXPECT_EQ(exitCode, EXIT_SUCCESS); +} + +TEST_F(LauncherTest, launches_multiObject_singleObject) +{ + int ac = 3; + const char* av[] = { NULL, "-i", "testData/cases/multiObject/singleObject.tessellator.json"}; + int exitCode; + EXPECT_NO_THROW(exitCode = launcher(ac, av)); + EXPECT_EQ(exitCode, EXIT_SUCCESS); +} + +TEST_F(LauncherTest, launches_multiObject_sameFileMultipleGroups) +{ + int ac = 3; + const char* av[] = { NULL, "-i", "testData/cases/multiObject/sameFileMultipleGroups.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 712fbfa..4e89dc5 100644 --- a/test/app/vtkIOTest.cpp +++ b/test/app/vtkIOTest.cpp @@ -23,8 +23,8 @@ TEST_F(VTKIOTest, readMeshFromSTL) TEST_F(VTKIOTest, exportAndReadMeshFromVTU) { auto mSTL{ readInputMesh("testData/cases/alhambra/alhambra.stl") }; - exportMeshToVTU("tmp_exported_alhambra.vtu", mSTL); - auto mVTU{ readInputMesh("tmp_exported_alhambra.vtu") }; + exportMeshToVTU("testData/cases/alhambra/tmp_exported_alhambra.vtu", mSTL); + auto mVTU{ readInputMesh("testData/cases/alhambra/tmp_exported_alhambra.vtu") }; EXPECT_EQ(mSTL.coordinates.size(), mVTU.coordinates.size()); EXPECT_EQ(mSTL.groups.size(), mVTU.groups.size()); @@ -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/CollapserTest.cpp b/test/core/CollapserTest.cpp index 126f887..e81c0ec 100644 --- a/test/core/CollapserTest.cpp +++ b/test/core/CollapserTest.cpp @@ -7,7 +7,10 @@ #include "utils/CoordGraph.h" #include "utils/GridTools.h" #include "utils/MeshTools.h" -#include "app/vtkIO.h" + +#if APP_LOADED + #include "app/vtkIO.h" +#endif namespace meshlib::core { @@ -146,6 +149,8 @@ TEST_F(CollapserTest, preserves_closedness) } +#if APP_LOADED + TEST_F(CollapserTest, closedness_for_sphere) { auto m = vtkIO::readInputMesh("testData/cases/sphere/sphere.stl"); @@ -164,6 +169,8 @@ TEST_F(CollapserTest, closedness_for_sphere) } +#endif + TEST_F(CollapserTest, areas_are_below_threshold_issue) { Mesh m; diff --git a/test/core/CompressorTest.cpp b/test/core/CompressorTest.cpp new file mode 100644 index 0000000..4a7a4f8 --- /dev/null +++ b/test/core/CompressorTest.cpp @@ -0,0 +1,503 @@ +#include +#include "core/Compressor.h" +#include "MeshFixtures.h" +#include "utils/MeshTools.h" + +using namespace meshlib; +using namespace meshlib::utils::meshTools; + +namespace meshlib::tests { + +class CompressorTest : public ::testing::Test { +protected: + void SetUp() override { + grid_ = { + std::vector{0, 1, 2, 3, 4, 5, 6}, + std::vector{0, 1, 2, 3, 4, 5, 6}, + std::vector{0, 1, 2, 3, 4, 5, 6} + }; + } + + Grid grid_; +}; + +TEST_F(CompressorTest, Compress2x2QuadsIntoOneSurface) { + // Create 4 quads arranged in a 2x2 pattern on the same plane + // They should be merged into a single surface + + Mesh mesh; + mesh.grid = grid_; + + // Quad 1: bottom-left (cells 0,0 to 1,1) + addQuad(mesh, {0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}); + // Quad 2: bottom-right (cells 1,0 to 2,1) + addQuad(mesh, {1, 0, 0}, {2, 0, 0}, {2, 1, 0}, {1, 1, 0}); + // Quad 3: top-left (cells 0,1 to 1,2) + addQuad(mesh, {0, 1, 0}, {1, 1, 0}, {1, 2, 0}, {0, 2, 0}); + // Quad 4: top-right (cells 1,1 to 2,2) + addQuad(mesh, {1, 1, 0}, {2, 1, 0}, {2, 2, 0}, {1, 2, 0}); + + EXPECT_EQ(countMeshElementsIf(mesh, isQuad), 4u); + + auto merged = core::Compressor::compressSurfacesInMesh(mesh); + + EXPECT_EQ(merged, 3u); + ASSERT_EQ(mesh.groups.size(), 1u); + ASSERT_EQ(mesh.groups[0].elements.size(), 1u); + + EXPECT_EQ( + CoordinateIds({0, 4, 8, 7}), + mesh.groups[0].elements[0].vertices + ); + +} + +TEST_F(CompressorTest, CompresRepeatedQuadsIntoOneSurface) { + // Create 4 quads arranged in a 2x2 pattern on the same plane + // They should be merged into a single surface + + Mesh mesh; + mesh.grid = grid_; + + addQuad(mesh, {0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}); + addQuad(mesh, {0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}); + + EXPECT_EQ(countMeshElementsIf(mesh, isQuad), 2u); + + auto merged = core::Compressor::compressSurfacesInMesh(mesh); + + EXPECT_EQ(merged, 1u); + ASSERT_EQ(mesh.groups.size(), 1u); + ASSERT_EQ(mesh.groups[0].elements.size(), 1u); + + EXPECT_EQ( + CoordinateIds({0, 1, 2, 3}), + mesh.groups[0].elements[0].vertices + ); + +} + +TEST_F(CompressorTest, DoesNotCompressQuadsWithDifferentOrientation) { + // Create quads on different planes - should not be merged + + Mesh mesh; + mesh.grid = grid_; + + // Counter clock-wise quad + addQuad(mesh, {0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}); + // Adjacent clock-wise quad + addQuad(mesh, {1, 0, 0}, {1, 1, 0}, {2, 1, 0}, {1, 0, 0}); + + EXPECT_EQ(countMeshElementsIf(mesh, isQuad), 2u); + + auto merged = core::Compressor::compressSurfacesInMesh(mesh); + + EXPECT_EQ(merged, 0u); + EXPECT_EQ(countMeshElementsIf(mesh, isQuad), 2u); +} + +TEST_F(CompressorTest, DoesNotCompressNonCoplanarQuads) { + // Create quads on different planes - should not be merged + + Mesh mesh; + mesh.grid = grid_; + + // Quad on z=0 plane + addQuad(mesh, {0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}); + // Quad on z=1 plane (different plane) + addQuad(mesh, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}, {0, 1, 1}); + + EXPECT_EQ(countMeshElementsIf(mesh, isQuad), 2u); + + auto merged = core::Compressor::compressSurfacesInMesh(mesh); + + EXPECT_EQ(merged, 0u); + EXPECT_EQ(countMeshElementsIf(mesh, isQuad), 2u); +} + +TEST_F(CompressorTest, DoesNotCompressDisconnectedQuads) { + // Create quads on same plane but not connected - should not be merged + + Mesh mesh; + mesh.grid = grid_; + + // Quad at bottom-left + addQuad(mesh, {0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}); + // Quad at top-right (disconnected) + addQuad(mesh, {3, 3, 0}, {4, 3, 0}, {4, 4, 0}, {3, 4, 0}); + + EXPECT_EQ(countMeshElementsIf(mesh, isQuad), 2u); + + auto merged = core::Compressor::compressSurfacesInMesh(mesh); + + EXPECT_EQ(merged, 0u); + EXPECT_EQ(countMeshElementsIf(mesh, isQuad), 2u); +} + +TEST_F(CompressorTest, CompressWithHoleCreatesInnerContour) { + // Create 8 quads forming a ring with a hole in the middle + // The ring decomposes into 4 rectangles (left col, right col, top center, bottom center) + + Mesh mesh; + mesh.grid = grid_; + + // Outer ring of quads (leaving center 1,1 to 2,2 empty) + // Layout: + // x=0 x=1 x=2 x=3 + // y=3 + // [5] [8] [6] + // y=2 + // [3] hole [4] + // y=1 + // [1] [7] [2] + // y=0 + // Bottom row + addQuad(mesh, {0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}); // 0, 1, 2, 3 + addQuad(mesh, {2, 0, 0}, {3, 0, 0}, {3, 1, 0}, {2, 1, 0}); // 4, 5, 6, 7 + // Middle row (sides only) + addQuad(mesh, {0, 1, 0}, {1, 1, 0}, {1, 2, 0}, {0, 2, 0}); // 3, 2, 8, 9 + addQuad(mesh, {2, 1, 0}, {3, 1, 0}, {3, 2, 0}, {2, 2, 0}); // 7, 6, 10, 11 + // Top row + addQuad(mesh, {0, 2, 0}, {1, 2, 0}, {1, 3, 0}, {0, 3, 0}); // 9, 8, 12, 13 + addQuad(mesh, {2, 2, 0}, {3, 2, 0}, {3, 3, 0}, {2, 3, 0}); // 11, 10, 14, 15 + // Corners to complete the ring + addQuad(mesh, {1, 0, 0}, {2, 0, 0}, {2, 1, 0}, {1, 1, 0}); // 1, 4, 7, 2 + addQuad(mesh, {1, 3, 0}, {1, 2, 0}, {2, 2, 0}, {2, 3, 0}); // 12, 8, 11, 15 + + EXPECT_EQ(countMeshElementsIf(mesh, isQuad), 8u); + + auto merged = core::Compressor::compressSurfacesInMesh(mesh); + + auto finalCount = countMeshElementsIf(mesh, isQuad); + + // Optimal decomposition: 4 rectangles + // - Bottom row (quads 1,7,2): cells x=0-3, y=0-1 + // - Left center (quad 3): cell x=0-1, y=1-2 + // - Top row (quads 5, 8, 6): cell x=0-3, y=2-3 + // - Right center (quad 4): cell x=2-3, y=1-2 + EXPECT_EQ(finalCount, 4u); + EXPECT_EQ(merged, 4u); // 8 - 4 = 4 surfaces merged + + EXPECT_EQ( + CoordinateIds({0, 5, 6, 3}), + mesh.groups[0].elements[0].vertices + ); + + EXPECT_EQ( + CoordinateIds({3, 2, 8, 9}), + mesh.groups[0].elements[1].vertices + ); + + EXPECT_EQ( + CoordinateIds({9, 10, 14, 13}), + mesh.groups[0].elements[2].vertices + ); + + EXPECT_EQ( + CoordinateIds({7, 6, 10, 11}), + mesh.groups[0].elements[3].vertices + ); +} + +TEST_F(CompressorTest, DoesNotCompressQuadsWithDifferentNormals) { + // Create quads with different normal directions - should not be merged + + Mesh mesh; + mesh.grid = grid_; + + // Quad on z=0 plane, normal pointing +z + addQuad(mesh, {0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}); + // Quad on y=0 plane, normal pointing +y (different orientation) + addQuad(mesh, {0, 0, 0}, {1, 0, 0}, {1, 0, 1}, {0, 0, 1}); + + EXPECT_EQ(countMeshElementsIf(mesh, isQuad), 2u); + + auto merged = core::Compressor::compressSurfacesInMesh(mesh); + + EXPECT_EQ(merged, 0u); + EXPECT_EQ(countMeshElementsIf(mesh, isQuad), 2u); +} + +TEST_F(CompressorTest, CompressRingRoundTrip) { + // Create ring of 8 quads, compress + + Mesh mesh; + mesh.grid = grid_; + + // Outer ring of quads (leaving center 1,1 to 2,2 empty) + // Layout: + // x=0 x=1 x=2 x=3 x=4 + // y=4 + // [10] [09] [08] [07] + // y=3 + // [11] [06] + // y=2 + // [12] [05] + // y=1 + // [01] [02] [03] [04] + // y=0 + // Bottom row + + addQuad(mesh, {0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}); // 01 + addQuad(mesh, {1, 0, 0}, {2, 0, 0}, {2, 1, 0}, {1, 1, 0}); // 02 + addQuad(mesh, {2, 0, 0}, {3, 0, 0}, {3, 1, 0}, {2, 1, 0}); // 03 + addQuad(mesh, {3, 0, 0}, {4, 0, 0}, {4, 1, 0}, {3, 1, 0}); // 04 + addQuad(mesh, {3, 1, 0}, {4, 1, 0}, {4, 2, 0}, {3, 2, 0}); // 05 + addQuad(mesh, {3, 2, 0}, {4, 2, 0}, {4, 3, 0}, {3, 3, 0}); // 06 + addQuad(mesh, {3, 3, 0}, {4, 3, 0}, {4, 4, 0}, {3, 4, 0}); // 07 + addQuad(mesh, {2, 3, 0}, {3, 3, 0}, {3, 4, 0}, {2, 4, 0}); // 08 + addQuad(mesh, {1, 4, 0}, {1, 3, 0}, {2, 3, 0}, {2, 4, 0}); // 09 + addQuad(mesh, {0, 3, 0}, {1, 3, 0}, {1, 4, 0}, {0, 4, 0}); // 10 + addQuad(mesh, {0, 2, 0}, {1, 2, 0}, {1, 3, 0}, {0, 3, 0}); // 11 + addQuad(mesh, {0, 1, 0}, {1, 1, 0}, {1, 2, 0}, {0, 2, 0}); // 12 + + EXPECT_EQ(countMeshElementsIf(mesh, isQuad), 12u); + + // Compress: 12 quads -> 4 surfaces () + auto merged = core::Compressor::compressSurfacesInMesh(mesh); + EXPECT_EQ(merged, 8u); // 12 - 4 = 4 surfaces merged + auto compressedCount = countMeshElementsIf(mesh, isQuad); + EXPECT_EQ(compressedCount, 4u); +} + +TEST_F(CompressorTest, Compress3x3GridRoundTrip) { + // Create 9 quads in 3x3 grid, compress to 1 surface + + Mesh mesh; + mesh.grid = grid_; + + // 3x3 grid of quads + // Layout: + // x=0 x=1 x=2 x=3 + // y=3 + // [3] [6] [9] + // y=2 + // [2] [5] [8] + // y=1 + // [1] [4] [7] + // y=0 + // + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + addQuad(mesh, + {i, j, 0}, {i+1, j, 0}, + {i+1, j+1, 0}, {i, j+1, 0}); + } + } + + EXPECT_EQ(countMeshElementsIf(mesh, isQuad), 9u); + + // Compress: 9 quads -> 1 surface + auto merged = core::Compressor::compressSurfacesInMesh(mesh); + EXPECT_EQ(merged, 8u); + EXPECT_EQ(countMeshElementsIf(mesh, isQuad), 1u); + + EXPECT_EQ( + CoordinateIds({0, 12, 15, 7}), + mesh.groups[0].elements[0].vertices + ); +} + +// ============== Line Compression Tests ============== + +TEST_F(CompressorTest, Compress2CollinearLinesIntoOne) { + Mesh mesh; + mesh.grid = grid_; + + addLine(mesh, {0, 0, 0}, {1, 0, 0}); + addLine(mesh, {1, 0, 0}, {2, 0, 0}); + + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 2u); + + auto merged = core::Compressor::compressLinesInMesh(mesh); + + EXPECT_EQ(merged, 1u); + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 1u); + + ASSERT_EQ(mesh.groups[0].elements.size(), 1u); + const auto& line = mesh.groups[0].elements[0]; + EXPECT_EQ(line.type, Element::Type::Line); + EXPECT_EQ(CoordinateIds({0, 2}), line.vertices); +} + +TEST_F(CompressorTest, CompressRepeatedLinesIntoOne) { + Mesh mesh; + mesh.grid = grid_; + + addLine(mesh, {0, 0, 0}, {1, 0, 0}); + addLine(mesh, {0, 0, 0}, {1, 0, 0}); + + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 2u); + + auto merged = core::Compressor::compressLinesInMesh(mesh); + + EXPECT_EQ(merged, 1u); + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 1u); + + ASSERT_EQ(mesh.groups[0].elements.size(), 1u); + const auto& line = mesh.groups[0].elements[0]; + EXPECT_EQ(line.type, Element::Type::Line); + EXPECT_EQ(CoordinateIds({0, 1}), line.vertices); +} + +TEST_F(CompressorTest, DoesNotCompressNonCollinearLines) { + Mesh mesh; + mesh.grid = grid_; + + addLine(mesh, {0, 0, 0}, {1, 0, 0}); + addLine(mesh, {0, 0, 0}, {0, 1, 0}); + + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 2u); + + auto merged = core::Compressor::compressLinesInMesh(mesh); + + EXPECT_EQ(merged, 0u); + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 2u); +} + +TEST_F(CompressorTest, DoesNotCompressDisconnectedLines) { + Mesh mesh; + mesh.grid = grid_; + + addLine(mesh, {0, 0, 0}, {1, 0, 0}); + addLine(mesh, {3, 0, 0}, {4, 0, 0}); + + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 2u); + + auto merged = core::Compressor::compressLinesInMesh(mesh); + + EXPECT_EQ(merged, 0u); + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 2u); +} + +TEST_F(CompressorTest, DoesNotCompressOverlappingOppositeDirectionLines) { + Mesh mesh; + mesh.grid = grid_; + + addLine(mesh, {0, 0, 0}, {1, 0, 0}); + addLine(mesh, {1, 0, 0}, {0, 0, 0}); + + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 2u); + + auto merged = core::Compressor::compressLinesInMesh(mesh); + + EXPECT_EQ(merged, 0u); + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 2u); +} + +TEST_F(CompressorTest, DoesNotCompressConnectedOppositeDirectionLines) { + Mesh mesh; + mesh.grid = grid_; + + addLine(mesh, {0, 0, 0}, {1, 0, 0}); + addLine(mesh, {2, 0, 0}, {1, 0, 0}); + + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 2u); + + auto merged = core::Compressor::compressLinesInMesh(mesh); + + EXPECT_EQ(merged, 0u); + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 2u); +} + +TEST_F(CompressorTest, Compress3LinesIntoOne) { + Mesh mesh; + mesh.grid = grid_; + + addLine(mesh, {0, 0, 0}, {1, 0, 0}); + addLine(mesh, {1, 0, 0}, {2, 0, 0}); + addLine(mesh, {2, 0, 0}, {3, 0, 0}); + + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 3u); + + auto merged = core::Compressor::compressLinesInMesh(mesh); + + EXPECT_EQ(merged, 2u); + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 1u); + EXPECT_EQ(CoordinateIds({0, 3}), mesh.groups[0].elements[0].vertices); +} + +TEST_F(CompressorTest, Compress5LineRoundTrip) { + Mesh mesh; + mesh.grid = grid_; + + addLine(mesh, {0, 0, 0}, {0, 1, 0}); + addLine(mesh, {0, 1, 0}, {0, 2, 0}); + addLine(mesh, {0, 2, 0}, {0, 3, 0}); + addLine(mesh, {0, 3, 0}, {0, 4, 0}); + addLine(mesh, {0, 4, 0}, {0, 5, 0}); + + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 5u); + + core::Compressor::compressLinesInMesh(mesh); + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 1u); + EXPECT_EQ(CoordinateIds({0, 5}), mesh.groups[0].elements[0].vertices); +} + +TEST_F(CompressorTest, CompressMixedDirections) { + Mesh mesh; + mesh.grid = grid_; + + // Layout: + // z + // v6 + // |l6 + // v5 + // |l5 + // v0———v1———v2 x + // l3 / l1 l2 + // v3 + // l4 / + // v4 + // y + + addLine(mesh, {0, 0, 0}, {1, 0, 0}); + addLine(mesh, {1, 0, 0}, {2, 0, 0}); + addLine(mesh, {0, 0, 0}, {0, 1, 0}); + addLine(mesh, {0, 1, 0}, {0, 2, 0}); + addLine(mesh, {0, 0, 0}, {0, 0, 1}); + addLine(mesh, {0, 0, 1}, {0, 0, 2}); + + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 6u); + + auto merged = core::Compressor::compressLinesInMesh(mesh); + + EXPECT_EQ(merged, 3u); + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 3u); + + EXPECT_EQ(CoordinateIds({0, 2}), mesh.groups[0].elements[0].vertices); + EXPECT_EQ(CoordinateIds({0, 4}), mesh.groups[0].elements[1].vertices); + EXPECT_EQ(CoordinateIds({0, 6}), mesh.groups[0].elements[2].vertices); +} + +TEST_F(CompressorTest, DoesNotCompressLinesWhenGroupDimensionIsLowerThanSurface) { + Mesh mesh; + mesh.grid = grid_; + + mesh.groups.resize(3); + + for (GroupId g = 0; g < 3; ++g){ + addLine(mesh, {0, 0, 0}, {0, 1, 0}, g); + addLine(mesh, {0, 1, 0}, {0, 2, 0}, g); + addLine(mesh, {0, 2, 0}, {0, 3, 0}, g); + addLine(mesh, {0, 3, 0}, {0, 4, 0}, g); + addLine(mesh, {0, 4, 0}, {0, 5, 0}, g); + } + + auto dimensions = std::vector({Element::Type::Volume, Element::Type::Surface, Element::Type::Line}); + + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 15u); + + auto merged = core::Compressor::compressLinesInMesh(mesh, dimensions); + + EXPECT_EQ(merged, 8u); + EXPECT_EQ(countMeshElementsIf(mesh, isLine), 7u); + + EXPECT_EQ(1, mesh.groups[0].elements.size()); + EXPECT_EQ(1, mesh.groups[1].elements.size()); + EXPECT_EQ(5, mesh.groups[2].elements.size()); + EXPECT_EQ(CoordinateIds({0, 5}), mesh.groups[0].elements[0].vertices); + EXPECT_EQ(CoordinateIds({0, 5}), mesh.groups[1].elements[0].vertices); +} + +} diff --git a/test/core/SlicerTest.cpp b/test/core/SlicerTest.cpp index 1ae8c6e..2367b93 100644 --- a/test/core/SlicerTest.cpp +++ b/test/core/SlicerTest.cpp @@ -5,10 +5,13 @@ #include "Collapser.h" #include "Geometry.h" #include "MeshTools.h" -#include "app/vtkIO.h" #include "utils/RedundancyCleaner.h" #include "utils/CoordGraph.h" +#if APP_LOADED + #include "app/vtkIO.h" +#endif + namespace meshlib::core { using namespace meshFixtures; @@ -990,6 +993,7 @@ TEST_F(SlicerTest, canKeepOppositeLinesInLineDimensionPolicy) } } +#if APP_LOADED TEST_F(SlicerTest, preserves_topological_closedness_for_alhambra) { @@ -1024,6 +1028,8 @@ TEST_F(SlicerTest, preserves_topological_closedness_for_sphere) // vtkIO::exportMeshToVTU("testData/cases/sphere/sphere.contour.vtk", contourMesh); } +#endif + TEST_F(SlicerTest, sphere_case_patch_contour_check_1) { Mesh m; @@ -1044,14 +1050,17 @@ TEST_F(SlicerTest, sphere_case_patch_contour_check_1) EXPECT_EQ(countContours(m), countContours(slicedMesh)); + #if APP_LOADED //For debugging. // meshTools::convertToAbsoluteCoordinates(slicedMesh); // vtkIO::exportMeshToVTU("testData/cases/sphere/sphere.sliced.vtk", slicedMesh); // auto contourMesh = meshTools::buildMeshFromContours(slicedMesh); // vtkIO::exportMeshToVTU("testData/cases/sphere/sphere.contour.vtk", contourMesh); + #endif } + TEST_F(SlicerTest, sphere_case_patch_contour_check_2) { Mesh m; @@ -1071,13 +1080,14 @@ TEST_F(SlicerTest, sphere_case_patch_contour_check_2) auto slicedMesh = Slicer{m}.getMesh(); EXPECT_EQ(countContours(m), countContours(slicedMesh)); - + #if APP_LOADED // For debugging. - meshTools::convertToAbsoluteCoordinates(slicedMesh); - vtkIO::exportMeshToVTU("sliced.vtk", slicedMesh); + // meshTools::convertToAbsoluteCoordinates(slicedMesh); + // vtkIO::exportMeshToVTU("sliced.vtk", slicedMesh); - auto contourMesh = meshTools::buildMeshFromContours(slicedMesh); - vtkIO::exportMeshToVTU("contour.vtk", contourMesh); + // auto contourMesh = meshTools::buildMeshFromContours(slicedMesh); + // vtkIO::exportMeshToVTU("contour.vtk", contourMesh); + #endif } -} \ No newline at end of file +} diff --git a/test/core/SmootherTest.cpp b/test/core/SmootherTest.cpp index 228ddb0..fce3dac 100644 --- a/test/core/SmootherTest.cpp +++ b/test/core/SmootherTest.cpp @@ -6,7 +6,10 @@ #include "utils/Geometry.h" #include "utils/MeshTools.h" #include "core/Slicer.h" -#include "app/vtkIO.h" + +#if APP_LOADED + #include "app/vtkIO.h" +#endif namespace meshlib::core { using namespace utils; @@ -87,6 +90,8 @@ TEST_F(SmootherTest, touching_by_single_point) EXPECT_EQ(1, countMeshElementsIf(r, isTriangle)); } +#if APP_LOADED + TEST_F(SmootherTest, preserves_topological_closedness_for_alhambra) { @@ -116,7 +121,6 @@ TEST_F(SmootherTest, preserves_topological_closedness_for_alhambra) // vtkIO::exportMeshToVTU("testData/cases/alhambra/alhambra.contour.vtk", contourMesh); } - TEST_F(SmootherTest, preserves_topological_closedness_for_sphere) { auto m = vtkIO::readInputMesh("testData/cases/sphere/sphere.stl"); @@ -146,4 +150,6 @@ TEST_F(SmootherTest, preserves_topological_closedness_for_sphere) // vtkIO::exportMeshToVTU("testData/cases/sphere/sphere.contour.vtk", contourMesh); } +#endif + } \ No newline at end of file diff --git a/test/core/SmootherToolsTest.cpp b/test/core/SmootherToolsTest.cpp index 78e2670..e93a9d0 100644 --- a/test/core/SmootherToolsTest.cpp +++ b/test/core/SmootherToolsTest.cpp @@ -419,6 +419,30 @@ class SmootherToolsTest : public ::testing::Test { const double alignmentAngle = 5.0; }; +class SmootherToolsTestAccess { +public: + static CoordGraph::Path pathFromIdToAnyTarget( + const CoordinateId startId, + const CoordGraph::Path& cycle, + const bool forward, + const IdSet& target) + { + return SmootherTools::pathFromIdToAnyTarget(startId, cycle, forward, target); + } +}; + +TEST_F(SmootherToolsTest, pathFromLargeCoordinateIdContainsTheIdOnce) +{ + const CoordinateId startId = 1000; + const CoordinateId targetId = 1001; + const CoordGraph::Path cycle = {999, startId, targetId}; + + const auto path = SmootherToolsTestAccess::pathFromIdToAnyTarget( + startId, cycle, true, {targetId}); + + EXPECT_EQ((CoordGraph::Path{startId, targetId}), path); +} + /// Remesh patch with interior points to remove them. /// \verbatim /// 1 ---- 2 @@ -673,6 +697,105 @@ TEST_F(SmootherToolsTest, collapsePointsOnFeatureEdges_singlePatch) EXPECT_EQ(7, countDifferentCoordinates(cs)); } +TEST_F(SmootherToolsTest, collapsePointsOnFeatureEdges_withDisconnectedTargets) +{ + Mesh m; + m.grid = buildUnitLengthGrid(1.0); + m.coordinates = { + Coordinate({0.1, 0.1, 0.2}), + Coordinate({0.5, 0.1, 0.2}), + Coordinate({0.9, 0.1, 0.2}), + Coordinate({0.5, 0.4, 0.2}), + Coordinate({0.1, 0.6, 0.8}), + Coordinate({0.5, 0.6, 0.8}), + Coordinate({0.9, 0.6, 0.8}), + Coordinate({0.5, 0.9, 0.8}), + }; + m.groups = {Group()}; + m.groups[0].elements = { + Element({0, 1, 3}), + Element({1, 2, 3}), + Element({4, 5, 7}), + Element({5, 6, 7}), + }; + + Elements& elements = m.groups[0].elements; + const ElementsView disconnectedPatch = { + &elements[0], &elements[1], &elements[2], &elements[3] + }; + const SmootherTools::SingularIds singularIds( + {0, 1, 2, 4, 5, 6}, {}, {}); + + EXPECT_NO_THROW(SmootherTools(m.grid).collapsePointsOnFeatureEdges( + m.coordinates, disconnectedPatch, singularIds)); + EXPECT_EQ(m.coordinates[0], m.coordinates[1]); + EXPECT_EQ(m.coordinates[4], m.coordinates[5]); +} + +TEST_F(SmootherToolsTest, collapsePointsOnFeatureEdgesDoesNotCrossTouchingCell) +{ + Mesh m; + m.grid = buildUnitLengthGrid(0.25); + m.coordinates = { + Coordinate({0.9, 0.1, 0.2}), + Coordinate({1.2, 0.1, 0.2}), + Coordinate({1.9, 0.1, 0.2}), + Coordinate({1.5, 0.4, 0.2}), + }; + m.groups = {Group()}; + m.groups[0].elements = { + Element({0, 1, 3}), + Element({1, 2, 3}), + }; + + Elements& elements = m.groups[0].elements; + const ElementsView patch = {&elements[0], &elements[1]}; + const SmootherTools::SingularIds singularIds({0, 1, 2}, {}, {}); + + SmootherTools(m.grid).collapsePointsOnFeatureEdges( + m.coordinates, patch, singularIds); + + EXPECT_EQ(Coordinate({1.2, 0.1, 0.2}), m.coordinates[1]); + EXPECT_FALSE(GridTools(m.grid).elementCrossesGrid(elements[1], m.coordinates)); +} + +TEST_F(SmootherToolsTest, collapsePointsOnFeatureEdgesRejectsCombinedMovesThatCrossGrid) +{ + Mesh m; + m.grid = buildUnitLengthGrid(0.25); + m.coordinates = { + Coordinate({1.0, 0.2, 0.2}), + Coordinate({1.0, 1.0, 0.2}), + Coordinate({1.0, 1.0, 0.8}), + Coordinate({1.0, 1.8, 0.8}), + Coordinate({0.8, 1.0, 0.5}), + }; + m.groups = {Group()}; + m.groups[0].elements = { + Element({0, 1, 4}), + Element({1, 2, 4}), + Element({2, 3, 4}), + }; + + Elements& elements = m.groups[0].elements; + const ElementsView patch = {&elements[0], &elements[1], &elements[2]}; + const SmootherTools::SingularIds singularIds({0, 1, 2, 3}, {}, {}); + SmootherTools::IncidentElements incidentElements; + for (const auto& element : elements) { + for (const auto id : element.vertices) { + incidentElements[id].push_back(&element); + } + } + + SmootherTools(m.grid).collapsePointsOnFeatureEdges( + m.coordinates, patch, singularIds, incidentElements); + + EXPECT_EQ(m.coordinates[0], m.coordinates[1]); + EXPECT_NE(m.coordinates[3], m.coordinates[2]); + EXPECT_FALSE(GridTools(m.grid).elementCrossesGrid( + elements[1], m.coordinates)); +} + TEST_F(SmootherToolsTest, collapsePointsOnFeatureEdges_threePatches) { Mesh m = buildCornerMesh(); @@ -860,4 +983,4 @@ TEST_F(SmootherToolsTest, collapsePointsOnFeatureEdges_feature_in_interior) EXPECT_EQ(collapsed[5], collapsed[4]); } -} \ No newline at end of file +} diff --git a/test/core/SnapperTest.cpp b/test/core/SnapperTest.cpp index da42665..42cd629 100644 --- a/test/core/SnapperTest.cpp +++ b/test/core/SnapperTest.cpp @@ -1,13 +1,16 @@ #include "gtest/gtest.h" #include "MeshFixtures.h" - #include "Snapper.h" #include "Slicer.h" #include "Smoother.h" #include "utils/Tools.h" #include "utils/Geometry.h" #include "utils/MeshTools.h" -#include "app/vtkIO.h" + + +#if APP_LOADED + #include "app/vtkIO.h" +#endif using namespace meshlib; using namespace core; @@ -46,6 +49,8 @@ TEST_F(SnapperTest, similar_results_for_each_plane) } } +#if APP_LOADED + TEST_F(SnapperTest, preserves_topological_closedness_for_sphere) { auto m = vtkIO::readInputMesh("testData/cases/sphere/sphere.stl"); @@ -78,6 +83,7 @@ TEST_F(SnapperTest, preserves_topological_closedness_for_sphere) // vtkIO::exportMeshToVTU("testData/cases/sphere/sphere.contour.vtk", contourMesh); } +# endif // TEST_F(SnapperTest, triangles_convert_to_lines) // { diff --git a/test/core/StaircaserTest.cpp b/test/core/StaircaserTest.cpp index 4975e5f..baa63b1 100644 --- a/test/core/StaircaserTest.cpp +++ b/test/core/StaircaserTest.cpp @@ -2204,7 +2204,6 @@ TEST_F(StaircaserTest, selectiveStructurerWithEmptySetOfCells) } - TEST_F(StaircaserTest, modifyCoordinateOfASpecificCell) { 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/ConformalMesherTest.cpp b/test/meshers/ConformalMesherTest.cpp index bbc52a7..b13db5d 100644 --- a/test/meshers/ConformalMesherTest.cpp +++ b/test/meshers/ConformalMesherTest.cpp @@ -4,16 +4,24 @@ #include "meshers/ConformalMesher.h" #include "utils/Geometry.h" -#include "app/vtkIO.h" #include "utils/MeshTools.h" +#if APP_LOADED + #include "app/vtkIO.h" +#endif + namespace meshlib::meshers { using namespace meshFixtures; using namespace utils::meshTools; -using namespace vtkIO; + +#if APP_LOADED + using namespace vtkIO; +#endif class ConformalMesherTest : public ::testing::Test { protected: + +#if APP_LOADED Mesh launchConformalMesher(const std::string& inputFilename, const Mesh& inputMesh) { ConformalMesherOptions opts; @@ -31,6 +39,7 @@ class ConformalMesherTest : public ::testing::Test { return res; } +#endif }; TEST_F(ConformalMesherTest, cellsWithMoreThanAVertexPerEdge_1) @@ -333,6 +342,8 @@ TEST_F(ConformalMesherTest, cellsWithMoreThanAPathPerFace_8) EXPECT_EQ(2, res.size()); } +#if APP_LOADED + TEST_F(ConformalMesherTest, sphere) { // Input @@ -421,6 +432,8 @@ TEST_F(ConformalMesherTest, thinCylinder) EXPECT_NE(0, mesh.countElems()); } +#endif + // TEST_F(ConformalMesherTest, plane45_size05_grid_adapted) // { // ConformalMesher mesher(buildPlane45Mesh(0.5)); diff --git a/test/meshers/StaircaseMesherTest.cpp b/test/meshers/StaircaseMesherTest.cpp index 53d7373..2ec23d1 100644 --- a/test/meshers/StaircaseMesherTest.cpp +++ b/test/meshers/StaircaseMesherTest.cpp @@ -1,18 +1,26 @@ #include "gtest/gtest.h" #include "MeshFixtures.h" +#include +#include +#include +#include + #include "meshers/StaircaseMesher.h" +#include "StaircaseMesherOptions.h" #include "Staircaser.h" #include "core/Slicer.h" #include "core/Collapser.h" -#include "utils/Geometry.h" #include "utils/GridTools.h" #include "utils/MeshTools.h" #include "utils/RedundancyCleaner.h" -#include "app/vtkIO.h" + +#if APP_LOADED + #include "app/vtkIO.h" +#endif namespace meshlib::meshers { @@ -20,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: @@ -228,6 +362,7 @@ TEST_F(StaircaseMesherTest, testTriNonUniformGridStaircase) EXPECT_EQ(0, countMeshElementsIf(out, isNode)); } +#if APP_LOADED // FOR DEBUG ONLY / OBTAIN VISUAL REPRESENTATION TEST_F(StaircaseMesherTest, DISABLED_visualSelectiveStaircaserCone) @@ -265,16 +400,16 @@ TEST_F(StaircaseMesherTest, DISABLED_visualSelectiveStaircaserCone) } auto resultMesh = meshlib::core::Staircaser{ collapsedMesh }.getSelectiveMesh(cellSet); - // ASSERT_NO_THROW(meshTools::checkNoCellsAreCrossed(resultMesh)); + ASSERT_NO_THROW(meshTools::checkNoCellsAreCrossed(resultMesh)); RedundancyCleaner::removeOverlappedDimensionOneAndLowerElementsAndEquivalentSurfaces(resultMesh); utils::meshTools::reduceGrid(resultMesh, inputMesh.grid); utils::meshTools::convertToAbsoluteCoordinates(resultMesh); - // EXPECT_TRUE(meshTools::isAClosedTopology(inputMesh.groups[0].elements)); - // EXPECT_TRUE(meshTools::isAClosedTopology(surfaceMesh.groups[0].elements)); - // EXPECT_TRUE(meshTools::isAClosedTopology(slicedMesh.groups[0].elements)); - // EXPECT_TRUE(meshTools::isAClosedTopology(resultMesh.groups[0].elements)); + EXPECT_TRUE(meshTools::isAClosedTopology(inputMesh.groups[0].elements)); + EXPECT_TRUE(meshTools::isAClosedTopology(surfaceMesh.groups[0].elements)); + EXPECT_TRUE(meshTools::isAClosedTopology(slicedMesh.groups[0].elements)); + EXPECT_TRUE(meshTools::isAClosedTopology(resultMesh.groups[0].elements)); @@ -284,9 +419,9 @@ TEST_F(StaircaseMesherTest, DISABLED_visualSelectiveStaircaserCone) meshlib::vtkIO::exportGridToVTU(outputFolder / (basename + ".tessellator.selective.grid.vtk"), resultMesh.grid); } +#endif - -TEST_F(StaircaseMesherTest, DISABLED_testStaircaseTriangleWithUniformGrid) +TEST_F(StaircaseMesherTest, testStaircaseTriangleWithUniformGrid) { float lowerCoordinateValue = -0.5; @@ -311,10 +446,179 @@ TEST_F(StaircaseMesherTest, DISABLED_testStaircaseTriangleWithUniformGrid) ASSERT_NO_THROW(resultMesh = StaircaseMesher(inputMesh, 2).mesh()); EXPECT_EQ(0, countRepeatedElements(resultMesh)); - EXPECT_EQ(48, resultMesh.groups[0].elements.size()); + EXPECT_EQ(12, resultMesh.groups[0].elements.size()); EXPECT_EQ(10, countMeshElementsIf(resultMesh, isQuad)); - EXPECT_EQ(32, countMeshElementsIf(resultMesh, isLine)); - EXPECT_EQ(6, countMeshElementsIf(resultMesh, isNode)); + EXPECT_EQ(2, countMeshElementsIf(resultMesh, isLine)); + EXPECT_EQ(0, countMeshElementsIf(resultMesh, isNode)); +} + +TEST_F(StaircaseMesherTest, testStaircaseWithCompression) +{ + float lowerCoordinateValue = -0.5; + float upperCoordinateValue = 0.5; + int numberOfCells = 4; + float step = 0.25; + assert((upperCoordinateValue - lowerCoordinateValue) / (numberOfCells) == step); + + Mesh inputMesh; + inputMesh.grid = GridTools::buildCartesianGrid(lowerCoordinateValue, upperCoordinateValue, numberOfCells + 1); + inputMesh.coordinates = { + Coordinate({ -0.475, -0.15, 0 }), + Coordinate({ 0.475, -0.15, 0 }), + Coordinate({ 0.0 , 0.15, 0 }), + Coordinate({ 0.0 , 0.15, 0.475 }) + }; + inputMesh.groups.resize(2); + inputMesh.groups[0].elements = { + Element({0, 1, 2}, Element::Type::Surface), + Element({2, 3}, Element::Type::Line) + }; + inputMesh.groups[1].elements = { + Element({0, 2}, Element::Type::Line), + Element({2, 3}, Element::Type::Line) + }; + + Mesh nonCompressedMesh = StaircaseMesher(inputMesh, 2).mesh(); + Mesh compressedMesh; + StaircaseMesherOptions compressOption; + compressOption.compress = true; + ASSERT_NO_THROW(compressedMesh = StaircaseMesher(inputMesh, 2, compressOption).mesh()); + + EXPECT_EQ(3, countRepeatedElements(nonCompressedMesh)); + EXPECT_EQ(7, nonCompressedMesh.groups[0].elements.size()); + EXPECT_EQ(6, nonCompressedMesh.groups[1].elements.size()); + EXPECT_EQ(4, countMeshElementsIf(nonCompressedMesh, isQuad)); + EXPECT_EQ(9, countMeshElementsIf(nonCompressedMesh, isLine)); + EXPECT_EQ(0, countMeshElementsIf(nonCompressedMesh, isNode)); + + EXPECT_EQ(1, countRepeatedElements(compressedMesh)); + EXPECT_EQ(3, compressedMesh.groups[0].elements.size()); + EXPECT_EQ(6, nonCompressedMesh.groups[1].elements.size()); + EXPECT_EQ(1, countMeshElementsIf(compressedMesh, isQuad)); + EXPECT_EQ(8, countMeshElementsIf(compressedMesh, isLine)); + 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, fillsSphereAsSingleClosedUnitHexahedralVolume) +{ + auto mesh = vtkIO::readInputMesh("testData/cases/sphere/sphere.stl"); + + mesh.grid[X] = utils::GridTools::linspace(-100.0, 100.0, 51); + mesh.grid[Y] = utils::GridTools::linspace(-100.0, 100.0, 51); + mesh.grid[Z] = utils::GridTools::linspace(-100.0, 100.0, 51); + + ASSERT_TRUE(isSingleClosedSurface(mesh.groups[0].elements)); + + StaircaseMesherOptions options; + options.volumeGroups.insert(0); + options.splitHexahedra = true; + + // 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)); + } TEST_F(StaircaseMesherTest, preserves_topological_closedness_for_alhambra) @@ -499,5 +803,7 @@ TEST_F(StaircaseMesherTest, staircaser_reads_wires_correctly) // meshlib::vtkIO::exportGridToVTU(outputFolder / (basename + ".tessellator.selective.grid.vtk"), resultMesh.grid); } -} +#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; diff --git a/test/utils/CoordGraphTest.cpp b/test/utils/CoordGraphTest.cpp index 6a98958..54e9259 100644 --- a/test/utils/CoordGraphTest.cpp +++ b/test/utils/CoordGraphTest.cpp @@ -341,6 +341,20 @@ TEST_F(CoordGraphTest, ctors) EXPECT_EQ(cG.verticesSize(), cGViews.verticesSize()); EXPECT_EQ(cG.edgesSize(), cGViews.edgesSize()); } + +TEST_F(CoordGraphTest, constructorsAcceptSingleVertexInputs) +{ + const Elements elements = { + Element({7}, Element::Type::Node), + }; + const CoordGraph fromElements(elements); + const CoordGraph fromPaths(CoordGraph::Paths{{9}}); + + EXPECT_EQ(IdSet({7}), fromElements.getVertices()); + EXPECT_EQ(IdSet({9}), fromPaths.getVertices()); + EXPECT_EQ(0, fromElements.edgesSize()); + EXPECT_EQ(0, fromPaths.edgesSize()); +} TEST_F(CoordGraphTest, adding_edge_with_non_existingvertices) { CoordGraph g; @@ -563,6 +577,17 @@ TEST_F(CoordGraphTest, getClosestVerticesInSet_2) EXPECT_EQ(IdSet({ 3 }), g.getClosestVerticesInSet(2, { 3, 4 })); } + +TEST_F(CoordGraphTest, getClosestVerticesInSet_ignoresDisconnectedTargets) +{ + CoordGraph g; + g.addEdge(1, 2); + g.addEdge(3, 4); + + EXPECT_EQ(IdSet({2}), g.getClosestVerticesInSet(1, {2, 3})); + EXPECT_TRUE(g.getClosestVerticesInSet(1, {3}).empty()); +} + TEST_F(CoordGraphTest, graphIntersection) { @@ -1341,4 +1366,4 @@ TEST_F(CoordGraphTest, difference) } } -} \ No newline at end of file +} diff --git a/test/utils/MeshToolsTest.cpp b/test/utils/MeshToolsTest.cpp index b46c30a..8078d16 100644 --- a/test/utils/MeshToolsTest.cpp +++ b/test/utils/MeshToolsTest.cpp @@ -720,4 +720,90 @@ TEST_F(MeshToolsTest, reduceGrid_epsilon_coord) } } +TEST_F(MeshToolsTest, extractGroupsByName_singleGroup) +{ + Mesh m; + m.coordinates = { + Coordinate({0.0, 0.0, 0.0}), + Coordinate({1.0, 0.0, 0.0}), + Coordinate({0.0, 1.0, 0.0}) + }; + m.groups = { + Group("group_a", {Element({0, 1, 2}, Element::Type::Surface)}), + Group("group_b", {Element({0, 1}, Element::Type::Line)}) + }; + + Mesh result = extractGroupsByName(m, {"group_a"}); + + EXPECT_EQ(result.groups.size(), 1); + EXPECT_EQ(result.groups[0].name, "group_a"); + EXPECT_EQ(result.groups[0].elements.size(), 1); + EXPECT_EQ(result.coordinates.size(), 3); +} + +TEST_F(MeshToolsTest, extractGroupsByName_multipleGroups) +{ + Mesh m; + m.coordinates = { + Coordinate({0.0, 0.0, 0.0}), + Coordinate({1.0, 0.0, 0.0}), + Coordinate({0.0, 1.0, 0.0}), + Coordinate({1.0, 1.0, 0.0}) + }; + m.groups = { + Group("group_a", {Element({0, 1, 2}, Element::Type::Surface)}), + Group("group_b", {Element({0, 1}, Element::Type::Line)}), + Group("group_c", {Element({2, 3}, Element::Type::Line)}) + }; + + Mesh result = extractGroupsByName(m, {"group_a", "group_c"}); + + EXPECT_EQ(result.groups.size(), 2); + EXPECT_EQ(result.groups[0].name, "group_a"); + EXPECT_EQ(result.groups[0].elements.size(), 1); + EXPECT_EQ(result.groups[1].name, "group_c"); + EXPECT_EQ(result.groups[1].elements.size(), 1); +} + +TEST_F(MeshToolsTest, extractGroupsByName_nonExistentGroup) +{ + Mesh m; + m.coordinates = { + Coordinate({0.0, 0.0, 0.0}), + Coordinate({1.0, 0.0, 0.0}) + }; + m.groups = { + Group("group_a", {Element({0, 1}, Element::Type::Line)}) + }; + + Mesh result = extractGroupsByName(m, {"group_a", "non_existent"}); + + EXPECT_EQ(result.groups.size(), 2); + EXPECT_EQ(result.groups[0].name, "group_a"); + EXPECT_EQ(result.groups[0].elements.size(), 1); + EXPECT_EQ(result.groups[1].name, "non_existent"); + EXPECT_EQ(result.groups[1].elements.size(), 0); +} + +TEST_F(MeshToolsTest, extractGroupsByName_coordinateRemapping) +{ + Mesh m; + m.coordinates = { + Coordinate({0.0, 0.0, 0.0}), + Coordinate({1.0, 0.0, 0.0}), + Coordinate({0.0, 1.0, 0.0}) + }; + m.groups = { + Group("group_a", {Element({0, 1, 2}, Element::Type::Surface)}), + Group("group_b", {Element({0, 1}, Element::Type::Line)}) + }; + + Mesh result = extractGroupsByName(m, {"group_a", "group_b"}); + + EXPECT_EQ(result.groups.size(), 2); + EXPECT_EQ(result.groups[0].elements[0].vertices, std::vector({0, 1, 2})); + EXPECT_EQ(result.groups[1].elements[0].vertices, std::vector({3, 4})); + EXPECT_EQ(result.coordinates.size(), 5); +} + } \ No newline at end of file diff --git a/testData/cases/alhambra/alhambra.conformal.tessellator.json b/testData/cases/alhambra/alhambra.conformal.tessellator.json index f2b6c20..9f17e14 100644 --- a/testData/cases/alhambra/alhambra.conformal.tessellator.json +++ b/testData/cases/alhambra/alhambra.conformal.tessellator.json @@ -11,7 +11,8 @@ "type" : "conformal", "options" : { "edgePoints" : 6, - "forbiddenLength" : 0.15 + "forbiddenLength" : 0.15, + "compress": false } } } \ No newline at end of file diff --git a/testData/cases/alhambra/alhambra.tessellator.json b/testData/cases/alhambra/alhambra.tessellator.json index 2519ae2..f1d9032 100644 --- a/testData/cases/alhambra/alhambra.tessellator.json +++ b/testData/cases/alhambra/alhambra.tessellator.json @@ -6,5 +6,11 @@ [ 60.0, 60.0, 10.0] ] }, - "object": {"filename": "alhambra.stl"} + "object": {"filename": "alhambra.stl"}, + "mesher": { + "type" : "staircase", + "options" : { + "compress": true + } + } } \ No newline at end of file diff --git a/testData/cases/longPolyline/longPolyline.tessellator.json b/testData/cases/longPolyline/longPolyline.tessellator.json index 1c672f3..7f60504 100644 --- a/testData/cases/longPolyline/longPolyline.tessellator.json +++ b/testData/cases/longPolyline/longPolyline.tessellator.json @@ -6,5 +6,11 @@ [1100, 1100, 1100] ] }, - "object": {"filename": "longPolyline.vtu"} + "object": {"filename": "longPolyline.vtu"}, + "mesher": { + "type" : "staircase", + "options" : { + "compress" : false + } + } } \ No newline at end of file diff --git a/testData/cases/longPolyline/longPolyline_compression.tessellator.json b/testData/cases/longPolyline/longPolyline_compression.tessellator.json new file mode 100644 index 0000000..7294f08 --- /dev/null +++ b/testData/cases/longPolyline/longPolyline_compression.tessellator.json @@ -0,0 +1,16 @@ +{ + "grid": { + "numberOfCells": [20, 20, 20], + "boundingBox": [ + [ 600, 600, 600], + [1100, 1100, 1100] + ] + }, + "object": {"filename": "longPolyline.vtu"}, + "mesher": { + "type" : "staircase", + "options" : { + "compress" : true + } + } +} \ No newline at end of file diff --git a/testData/cases/longPolyline/longPolyline_legacy.tessellator.json b/testData/cases/longPolyline/longPolyline_legacy.tessellator.json new file mode 100644 index 0000000..1c672f3 --- /dev/null +++ b/testData/cases/longPolyline/longPolyline_legacy.tessellator.json @@ -0,0 +1,10 @@ +{ + "grid": { + "numberOfCells": [20, 20, 20], + "boundingBox": [ + [ 600, 600, 600], + [1100, 1100, 1100] + ] + }, + "object": {"filename": "longPolyline.vtu"} +} \ No newline at end of file diff --git a/testData/cases/multiObject/basic.tessellator.json b/testData/cases/multiObject/basic.tessellator.json new file mode 100644 index 0000000..5724c92 --- /dev/null +++ b/testData/cases/multiObject/basic.tessellator.json @@ -0,0 +1,16 @@ +{ + "grid": { + "numberOfCells": [50, 50, 50], + "boundingBox": [ + [-100.0, -100.0, -100.0], + [ 100.0, 100.0, 100.0] + ] + }, + "mesher": { + "type": "staircase" + }, + "objects": [ + {"filename": "sphere.stl", "volume": true, "group": "sphere_group"}, + {"filename": "cone.stl", "group": "cone_group"} + ] +} diff --git a/testData/cases/multiObject/cone.stl b/testData/cases/multiObject/cone.stl new file mode 100644 index 0000000..b855cac Binary files /dev/null and b/testData/cases/multiObject/cone.stl differ diff --git a/testData/cases/multiObject/mixedMesher.tessellator.json b/testData/cases/multiObject/mixedMesher.tessellator.json new file mode 100644 index 0000000..ca2afa7 --- /dev/null +++ b/testData/cases/multiObject/mixedMesher.tessellator.json @@ -0,0 +1,17 @@ +{ + "grid": { + "numberOfCells": [50, 50, 50], + "boundingBox": [ + [-100.0, -100.0, -100.0], + [ 100.0, 100.0, 100.0] + ] + }, + "mesher": { + "type": "staircase" + }, + "objects": [ + {"filename": "sphere.stl", "group": "sphere_group", "mesher": {"type": "staircase", "options": { "compress" : true} } }, + {"filename": "cone.stl", "group": "cone_group", "mesher": {"type": "conformal"}}, + {"filename": "cone.stl", "group": "default_group"} + ] +} diff --git a/testData/cases/multiObject/sameFileMultipleGroups.tessellator.json b/testData/cases/multiObject/sameFileMultipleGroups.tessellator.json new file mode 100644 index 0000000..b1bcd39 --- /dev/null +++ b/testData/cases/multiObject/sameFileMultipleGroups.tessellator.json @@ -0,0 +1,16 @@ +{ + "grid": { + "numberOfCells": [50, 50, 50], + "boundingBox": [ + [-100.0, -100.0, -100.0], + [ 100.0, 100.0, 100.0] + ] + }, + "mesher": { + "type": "staircase" + }, + "objects": [ + {"filename": "sphere.stl", "group": "group_a"}, + {"filename": "sphere.stl", "group": "group_b"} + ] +} diff --git a/testData/cases/multiObject/singleObject.tessellator.json b/testData/cases/multiObject/singleObject.tessellator.json new file mode 100644 index 0000000..9216505 --- /dev/null +++ b/testData/cases/multiObject/singleObject.tessellator.json @@ -0,0 +1,15 @@ +{ + "grid": { + "numberOfCells": [50, 50, 50], + "boundingBox": [ + [-100.0, -100.0, -100.0], + [ 100.0, 100.0, 100.0] + ] + }, + "mesher": { + "type": "staircase" + }, + "objects": [ + {"filename": "sphere.stl", "group": "default_group"} + ] +} diff --git a/testData/cases/multiObject/sphere.stl b/testData/cases/multiObject/sphere.stl new file mode 100644 index 0000000..29885a1 Binary files /dev/null and b/testData/cases/multiObject/sphere.stl differ diff --git a/testData/cases/sphere/closed_sphere.tessellator.json b/testData/cases/sphere/closed_sphere.tessellator.json new file mode 100644 index 0000000..d95ebfe --- /dev/null +++ b/testData/cases/sphere/closed_sphere.tessellator.json @@ -0,0 +1,10 @@ +{ + "grid": { + "numberOfCells": [50, 50, 50], + "boundingBox": [ + [-100.0, -100.0, -100.0], + [ 100.0, 100.0, 100.0] + ] + }, + "object": {"filename": "sphere.stl", "volume" : true} +} \ No newline at end of file diff --git a/vcpkg-configuration.json b/vcpkg-configuration.json index 563e537..e1d933f 100644 --- a/vcpkg-configuration.json +++ b/vcpkg-configuration.json @@ -10,5 +10,8 @@ "location": "https://aka.ms/vcpkg-ce-default", "name": "microsoft" } + ], + "overlay-ports" : [ + "./overlayPorts" ] }