diff --git a/.clang-format b/.clang-format new file mode 100644 index 0000000..db09510 --- /dev/null +++ b/.clang-format @@ -0,0 +1,32 @@ +--- +# Default rules +BasedOnStyle: LLVM +ColumnLimit: 120 +IndentWidth: 4 +UseTab: Never +--- +# C++ rules +Language: Cpp + +AccessModifierOffset: -4 +AlignAfterOpenBracket: DontAlign +AlignEscapedNewlines: Left +AllowAllArgumentsOnNextLine: false +AllowAllConstructorInitializersOnNextLine: false +AllowAllParametersOfDeclarationOnNextLine: false +AllowShortBlocksOnASingleLine: Empty +AllowShortFunctionsOnASingleLine: Empty +AllowShortIfStatementsOnASingleLine: Never +AllowShortLoopsOnASingleLine: false +AlwaysBreakTemplateDeclarations: Yes +BreakBeforeBraces: Attach +BreakConstructorInitializers: BeforeColon +ConstructorInitializerAllOnOneLineOrOnePerLine: true +DerivePointerAlignment: false +IndentCaseLabels: true +NamespaceIndentation: All +PointerAlignment: Left +SpaceAfterCStyleCast: true +SpaceAfterTemplateKeyword: false +Standard: c++17 +... diff --git a/.editorconfig b/.editorconfig new file mode 100644 index 0000000..409c268 --- /dev/null +++ b/.editorconfig @@ -0,0 +1,17 @@ +# EditorConfig: https://editorconfig.org + +root = true + +[*] +charset = utf-8 +indent_style = space +indent_size = 4 +insert_final_newline = true +trim_trailing_whitespace = true + +[{CMakeLists.txt,*.cmake}] +indent_size = 2 + +[*.md] +max_line_length = off +trim_trailing_whitespace = false diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml new file mode 100644 index 0000000..1dbd54f --- /dev/null +++ b/.github/workflows/build.yml @@ -0,0 +1,45 @@ +name: Build + +on: + push: + branches: [ main ] + +jobs: + linux: + name: Linux-Build + runs-on: ubuntu-22.04 + container: + image: kitware/paraview_org-plugin-devel:5.12.0 + env: + ACTIONS_ALLOW_USE_UNSECURE_NODE_VERSION: true # nodejs 16 is required for CentOS 7 + steps: + - name: Install dependencies + run: | + # Fix CentOS 7 EOL (https://serverfault.com/a/1161847) + sed -i s/mirror.centos.org/vault.centos.org/g /etc/yum.repos.d/CentOS-*.repo + sed -i s/^#.*baseurl=http/baseurl=http/g /etc/yum.repos.d/CentOS-*.repo + sed -i s/^mirrorlist=http/#mirrorlist=http/g /etc/yum.repos.d/CentOS-*.repo + yum -y install devtoolset-9 + - uses: actions/checkout@v3 # glibc in container is too old for v4 + - name: Configure + run: | + PARAVIEW_DIR=`find /builds/gitlab-kitware-sciviz-ci/build/install/lib/cmake -type d -regex "/builds/gitlab-kitware-sciviz-ci/build/install/lib/cmake/paraview-[0-9,.]*"` + scl enable devtoolset-9 -- cmake -S . -B ${{github.workspace}}/build -DParaView_DIR=$PARAVIEW_DIR -DCMAKE_BUILD_TYPE=Release + - name: Build + run: | + scl enable devtoolset-9 -- cmake --build ${{github.workspace}}/build --parallel 4 + - name: Install + run: | + scl enable devtoolset-9 -- cmake --install ${{github.workspace}}/build --prefix ./install + - name: Cleanup + run: | + mv ./install/lib64 ./install/lib + - uses: actions/upload-artifact@v3 # glibc in container is too old for v4 + with: + name: VofFlow + path: | + ./install/lib/paraview-5.12/plugins/VofFlow/VofFlow.so + ./install/bin/MaxVelocity + ./install/bin/PolyGap + ./install/bin/SeedGrid + ./install/bin/SmoothNormals diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..c692320 --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +build/ +cmake-build-*/ +.idea/ +.vs/ +.vscode/ diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..eabe4ee --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,81 @@ +cmake_minimum_required(VERSION 3.12...3.29 FATAL_ERROR) + +include(CMakeDependentOption) + +option(VOFFLOW_USE_VCPKG "Use vcpkg to load dependencies." ON) +cmake_dependent_option(VOFFLOW_DISABLE_VCPKG_BOOST "Disable vcpkg boost for ParaView superbuild." OFF "VOFFLOW_USE_VCPKG" OFF) + +# vcpkg +if (VOFFLOW_USE_VCPKG) + include(FetchContent) + mark_as_advanced(FORCE + FETCHCONTENT_BASE_DIR + FETCHCONTENT_FULLY_DISCONNECTED + FETCHCONTENT_QUIET + FETCHCONTENT_UPDATES_DISCONNECTED) + + # Require git for download + find_package(Git REQUIRED) + + FetchContent_Declare(vcpkg-download + GIT_REPOSITORY https://github.com/microsoft/vcpkg.git + GIT_TAG 2024.04.26 + SOURCE_DIR ${CMAKE_CURRENT_BINARY_DIR}/vcpkg) + FetchContent_GetProperties(vcpkg-download) + if (NOT vcpkg-download_POPULATED) + message(STATUS "Fetch vcpkg ...") + FetchContent_Populate(vcpkg-download) + mark_as_advanced(FORCE + FETCHCONTENT_SOURCE_DIR_VCPKG-DOWNLOAD + FETCHCONTENT_UPDATES_DISCONNECTED_VCPKG-DOWNLOAD) + endif () + + set(VCPKG_BOOTSTRAP_OPTIONS "-disableMetrics") + set(VCPKG_INSTALL_OPTIONS "--clean-after-build" "--no-print-usage") + set(CMAKE_TOOLCHAIN_FILE "${CMAKE_CURRENT_BINARY_DIR}/vcpkg/scripts/buildsystems/vcpkg.cmake" CACHE STRING "Vcpkg toolchain file") + set(ENV{VCPKG_FORCE_DOWNLOADED_BINARIES} ON) + + if (WIN32) + set(VCPKG_TARGET_TRIPLET "x64-windows-static") + endif () + + if (VOFFLOW_DISABLE_VCPKG_BOOST) + include(cmake/disable_vcpkg_boost.cmake) + endif () +endif () + +project(VofFlow + LANGUAGES C CXX) + +# Set a default build type if none was specified +if (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES) + message(STATUS "Setting build type to 'Release' as none was specified.") + set(CMAKE_BUILD_TYPE Release CACHE STRING "Choose the type of build." FORCE) + set_property(CACHE CMAKE_BUILD_TYPE PROPERTY STRINGS "Debug" "Release" "RelWithDebInfo" "MinSizeRel") +endif () + +find_package(ParaView REQUIRED) + +option(BUILD_SHARED_LIBS "Build shared libraries" ON) + +include(GNUInstallDirs) +set(CMAKE_RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_BINDIR}") +set(CMAKE_LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_LIBDIR}") +set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/${CMAKE_INSTALL_LIBDIR}") + +paraview_plugin_scan( + PLUGIN_FILES "${CMAKE_CURRENT_SOURCE_DIR}/src/ParaViewPlugin/paraview.plugin" + PROVIDES_PLUGINS plugins + ENABLE_BY_DEFAULT ON) + +paraview_plugin_build( + RUNTIME_DESTINATION "${CMAKE_INSTALL_BINDIR}" + LIBRARY_DESTINATION "${CMAKE_INSTALL_LIBDIR}" + LIBRARY_SUBDIRECTORY "${PARAVIEW_PLUGIN_SUBDIR}" + INSTALL_HEADERS OFF + PLUGINS ${plugins}) + +option(VOFFLOW_TOOLS "Build tools." ON) +if (VOFFLOW_TOOLS) + add_subdirectory(src/Tools) +endif () diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..df6e7a3 --- /dev/null +++ b/LICENSE @@ -0,0 +1,28 @@ +BSD 3-Clause License + +Copyright (c) 2024, Visualization Research Center (VISUS), University of Stuttgart + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + +1. Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. + +2. Redistributions in binary form must reproduce the above copyright notice, + this list of conditions and the following disclaimer in the documentation + and/or other materials provided with the distribution. + +3. Neither the name of the copyright holder nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE +FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL +DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR +SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER +CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, +OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/README.md b/README.md new file mode 100644 index 0000000..40a6423 --- /dev/null +++ b/README.md @@ -0,0 +1,120 @@ +# Visualization of Finite-Time Separation in Multiphase Flow + +*Moritz Heinemann, Johanna Potyka, Kathrin Schulte, Filip Sadlo, Thomas Ertl* + +![ParaView Plugin Properties](docs/img/teaser.png) + +This repository contains the ParaView plugin used in our [paper](https://doi.org/10.1109/TVCG.2024.3493607). +We provide [details on the reproducibility](docs/REPRODUCIBILITY.md) and usage examples in the [`docs`](docs) folder of this repository. + +## Download + +Plugin binaries for ParaView 5.12.0 are provided on the [Release Page](https://github.com/UniStuttgart-VISUS/vof-flow/releases). +They are compatible with the official ParaView 5.12.0 MPI release for Linux or Windows, available from the [ParaView website](https://www.paraview.org/download/): +- [Download ParaView 5.12.0 MPI Linux](https://www.paraview.org/paraview-downloads/download.php?submit=Download&version=v5.12&type=binary&os=Linux&downloadFile=ParaView-5.12.0-MPI-Linux-Python3.10-x86_64.tar.gz) +- [Download ParaView 5.12.0 MPI Windows](https://www.paraview.org/paraview-downloads/download.php?submit=Download&version=v5.12&type=binary&os=Windows&downloadFile=ParaView-5.12.0-MPI-Windows-Python3.10-msvc2017-AMD64.zip) + +## Build from source + +The plugin can be built from source using the common ParaView plugin interface. +For details, we refer to the [ParaView Plugin Howto](https://www.paraview.org/paraview-docs/v5.12.0/cxx/PluginHowto.html). +A more detailed example of the build steps can be found in our [reproducibility documentation](docs/REPRODUCIBILITY.md). + +## Example Data + +The jet-collision dataset from our paper can be used as example data: +[https://doi.org/10.18419/darus-4225](https://doi.org/10.18419/darus-4225). + +## Plugin Usage + +The plugin adds three filters to ParaView. +The main filter is called `VofTracking`, which implements the method described in our paper. +There are two additional filters called `Plic` and `Plic3` to generate the PLIC surfaces of Volume-of-Fluid datasets. +All filters require rectilinear grid data as input. +The properties of the `VofTracking` filters are documented below. + +### VofTracking Properties + +![ParaView Plugin Properties](docs/img/plugin_properties.png) + +| Property | Description | +|--------------------------|-------------------------------------------------------------------------------------------------------| +| `Use Three Phase` | Switch between two-phase and three-phase. | +| `VoF` | VoF array. | +| `VoF3` | VoF3 array (if three-phase). | +| `VoF Norm` | VoF normals array (if three-phase). | +| `Velocity` | Velocity array. | +| `Use Components` | Switch to use custom component fields. | +| `Components VoF` | Components field for the VoF data. | +| `Components VoF3` | Components field for the VoF3 data (if three-phase). | +| `Use Target Time Step` | Switch to use a fixed target time step. If false, the current ParaView time is used. | +| `Init Time Step` | Initial time step. | +| `Target Time Step` | Target time step (if `Use Target Time Step` is true, otherwise current ParaView time is used). | +| `Refinement` | Number of seeded particles per cell is (r+1)^3. | +| `Neighbor Correction` | Enable Neighbor Correction step. | +| `Cell Correction` | Enable Cell Correction step. | +| `PLIC Correction` | Enable PLIC Correction step. | +| `Integration Method`* | Select the particle advection integration method: `RK4` or `Euler`. | +| `Integration Sub-Steps`* | Select the number of sub-steps used for the integration. | +| `Epsilon`* | Epsilon for the PLIC surface generation. | +| `Num Iterations`* | Maximum number of iterations used for the PLIC surface generation. | +| `Ghost Cells` | Number of ghost cells shared between MPI processes. | +| `Cut Labels`* | Use label cut function (experimental). | +| `Label Cut Type`* | vtkImplicitFunction for cutting label space (experimental). | +| `Boundary Method`* | Method used to generate the separation boundaries `DiscreteMarchingCubes` or `DiscreteFlyingEdges3D`. | +| `Output Data Type` | Type of the grid output data: `vtkRectilinearGrid` or `vtkImageData`. | +| `Output State`* | Store state of plugin properties as json in the output data. | +| `Output Time Measure`* | Store benchmarking data in the output data. | +| `Mirror X Min` | Mirror boundary condition at X min boundary. | +| `Mirror X Max` | Mirror boundary condition at X max boundary. | +| `Mirror Y Min` | Mirror boundary condition at Y min boundary. | +| `Mirror Y Max` | Mirror boundary condition at Y max boundary. | +| `Mirror Z Min` | Mirror boundary condition at Z min boundary. | +| `Mirror Z Max` | Mirror boundary condition at Z max boundary. | + +Properties marked with `*` are advanced properties. +Use the ParaView `Toggle advanced properties` button to show them. + +### VofTracking Output + +![ParaView Plugin Properties](docs/img/plugin_pipeline.png) + +The VofTracking filter has four outputs. + +| Output | Description | +|--------------------|------------------------------------------------------------------------------------------------| +| Grid | Output grid with the extracted labels and the grid-based uncertainty values. | +| Seeds | Seeded particles, including label and particle-based uncertainty values. | +| Advected Particles | End position of the advected particles, including label and particle-based uncertainty values. | +| Boundaries | Extracted separation boundaries. | + +## Citing + +If you use our work in a scientific context, please cite our work: + +M. Heinemann, J. Potyka, K. Schulte, F. Sadlo, T. Ertl: "**Visualization of Finite-Time Separation in Multiphase Flow**" in IEEE Transactions on Visualization and Computer Graphics, doi: 10.1109/TVCG.2024.3493607. + +```bibtex +@Article{Heinemann2024VofFlow, + author = {Heinemann, Moritz and Potyka, Johanna and Schulte, Kathrin and Sadlo, Filip and Ertl, Thomas}, + journal = {IEEE Transactions on Visualization and Computer Graphics}, + title = {Visualization of Finite-Time Separation in Multiphase Flow}, + year = {2024}, + pages = {1-14}, + doi = {10.1109/TVCG.2024.3493607}, +} +``` + +## License + +> Copyright (c) 2024, Visualization Research Center (VISUS), University of Stuttgart + +Our code in this repository is licensed under the [BSD 3-Clause License](LICENSE). + +Additionally, the following dependencies with their corresponding licenses are used and/or downloaded automatically during the build step or in some of our scripts: + +- [ParaView](https://gitlab.kitware.com/paraview/paraview) +- [CGAL](https://github.com/cgal/cgal) +- [JSON for Modern C++](https://github.com/nlohmann/json) +- [vcpkg](https://github.com/microsoft/vcpkg) +- [ParaViewEasyPluginBuilder](https://gitlab.kitware.com/paraview/paraview-easy-plugin-builder) diff --git a/cmake/disable_vcpkg_boost.cmake b/cmake/disable_vcpkg_boost.cmake new file mode 100644 index 0000000..189cee5 --- /dev/null +++ b/cmake/disable_vcpkg_boost.cmake @@ -0,0 +1,21 @@ +# Disable all vcpkg boost* ports, to use the version provided by the ParaView superbuild. +# This overwrites all boost ports with empty ports. The empty ports use the original +# vcpkg.json to keep all features and dependencies next to an empty portfile. + +set(empty_ports_dir "${CMAKE_CURRENT_BINARY_DIR}/vcpkg_empty_ports") + +# Cleanup port dir +file(REMOVE_RECURSE "${empty_ports_dir}") +file(MAKE_DIRECTORY "${empty_ports_dir}") + +# Searching boost ports +set(port_dirs "${CMAKE_CURRENT_BINARY_DIR}/vcpkg/ports") +file(GLOB boost_ports RELATIVE ${port_dirs} ${port_dirs}/boost*) + +# Create empty port +foreach (portname ${boost_ports}) + file(COPY "${port_dirs}/${portname}/vcpkg.json" DESTINATION "${empty_ports_dir}/${portname}/") + file(WRITE "${empty_ports_dir}/${portname}/portfile.cmake" "set(VCPKG_POLICY_EMPTY_PACKAGE enabled)\n") +endforeach () + +set(VCPKG_OVERLAY_PORTS "${empty_ports_dir};${VCPKG_OVERLAY_PORTS}") diff --git a/docs/REPRODUCIBILITY.md b/docs/REPRODUCIBILITY.md new file mode 100644 index 0000000..b0ff2f8 --- /dev/null +++ b/docs/REPRODUCIBILITY.md @@ -0,0 +1,161 @@ +# Reproducibility Steps for Visualization of Finite-Time Separation in Multiphase Flow + +We tested this setup with Ubuntu 22.04. + +## Content + +- [Setup ParaView and the VofFlow Plugin](#setup-paraview-and-the-vofflow-plugin) +- [Download Dataset](#download-dataset) +- [Reproducing the Results](#reproducing-the-results) + +## Setup ParaView and the VofFlow Plugin + +### Option 1: Using Prebuild Binaries + +Plugin binaries for ParaView 5.12.0 are provided on the [Release Page](https://github.com/UniStuttgart-VISUS/vof-flow/releases). +They are compatible with the official ParaView 5.12.0 MPI release for Linux available from the [ParaView website](https://www.paraview.org/download/): + +- [Download ParaView 5.12.0 MPI Linux](https://www.paraview.org/paraview-downloads/download.php?submit=Download&version=v5.12&type=binary&os=Linux&downloadFile=ParaView-5.12.0-MPI-Linux-Python3.10-x86_64.tar.gz) + +Please unpack all downloaded files to `~/vofflow/install`, i.e., the `paraview` binary should be located at `~/vofflow/install/bin/paraview`. +This path is assumed in our example scripts. +Further, make the tool binaries from our plugin executable: `cd ~/vofflow/install/bin/ && chmod +x SeedGrid PolyGap SmoothNormals MaxVelocity`. + +> The released plugin binaries for Linux are compiled using the [GitHub Actions workflow in this repository](../.github/workflows/build.yml). +> The workflow is based on the [ParaViewEasyPluginBuilder](https://gitlab.kitware.com/paraview/paraview-easy-plugin-builder). +> The Windows version was compiled manually using the [ParaView Plugin Windows Binary Compatible Guide](https://gitlab.kitware.com/paraview/paraview-plugin-windows-binary-compatible-guide). + +### Option 2: Build Plugin using the ParaViewEasyPluginBuilder + +This option also uses the official ParaView release but compiles the plugin locally. +Please follow the steps from Option 1 to download and unpack the ParaView 5.12.0 release. +Please make sure Docker is installed on the system, see [https://docs.docker.com/engine/install/ubuntu/](https://docs.docker.com/engine/install/ubuntu/). + +Then, run [`./build_plugin_linux.sh`](scripts/build_plugin_linux.sh) to build the plugin binaries. +The plugin will be created in a folder called `plugin_build` in the current directory. +Please move the content from `plugin_build` to `~/vofflow/install`. + +> The script is based on the [ParaViewEasyPluginBuilder](https://gitlab.kitware.com/paraview/paraview-easy-plugin-builder). +> This is a Docker-based script to build ParaView plugins compatible with the official binary release of ParaView and does not require to build ParaView manually. +> We slightly modified the ParaViewEasyPluginBuilder by installing `devtoolset-9` for a C++17 compatible compiler. + +### Option 3: Compile ParaView and the Plugin (Ubuntu 22.04) + +- Install build environment: + ```shell + sudo apt install build-essential git cmake ninja-build + ``` +- Install ParaView dependencies: + ```shell + sudo apt install qtbase5-dev qttools5-dev libqt5svg5-dev qtxmlpatterns5-dev-tools libgl1-mesa-dev python3-dev libopenmpi-dev + ``` +- Install vcpkg dependencies: + ```shell + sudo apt install curl zip unzip tar pkg-config + ``` +- Clone ParaView and VofFlow: + ```shell + mkdir ~/vofflow && cd ~/vofflow + git clone --recursive --branch v5.12.0 https://gitlab.kitware.com/paraview/paraview.git + git clone https://github.com/UniStuttgart-VISUS/vof-flow.git + ``` +- Configure, Build, and Install ParaView: + ```shell + mkdir pv-build && cd pv-build + cmake -G Ninja -DCMAKE_BUILD_TYPE=Release -DPARAVIEW_USE_PYTHON=ON -DPARAVIEW_USE_MPI=ON -DCMAKE_INSTALL_PREFIX=../install ../paraview + cmake --build . + cmake --install . + cd .. + ``` +- Configure, Build, and Install VofFlow: + ```shell + mkdir vf-build && cd vf-build + cmake -G Ninja -DCMAKE_BUILD_TYPE=Release -DCMAKE_PREFIX_PATH=~/vofflow/install -DCMAKE_INSTALL_PREFIX=../install ../vof-flow/ + cmake --build . + cmake --install . + cd .. + ``` + +> This manual setup was tested with ParaView 5.11.1 and 5.12.0. + +#### Alternative using system packages instead of vcpkg + +By default, all dependencies except ParaView itself are downloaded automatically using vcpkg during the CMake configure of our project. +If you prefer system packages over vcpkg, you can install them manually: +```shell +sudo apt install libcgal-dev nlohmann-json3-dev +``` +In the VofFlow configure step, pass the additional parameter `-DVOFFLOW_USE_VCPKG=OFF`: +```shell +cmake -G Ninja -DCMAKE_BUILD_TYPE=Release -DCMAKE_PREFIX_PATH=~/vofflow/install -DCMAKE_INSTALL_PREFIX=../install -DVOFFLOW_USE_VCPKG=OFF ../vof-flow/ +``` + +#### Alternative for the ParaView-Superbuild + +If the ParaView-Superbuild is used to build ParaView, there could be version conflicts between the Boost version of the ParaView-Superbuild and the vcpkg provided Boost version. +We provide a CMake option to disable downloading Boost with vcpkg: `-DVOFFLOW_DISABLE_VCPKG_BOOST=ON`. + +## Download Dataset + +We have published the jet-collision dataset from our paper here: [https://doi.org/10.18419/darus-4225](https://doi.org/10.18419/darus-4225). +We provide the script `download_data.py` to download all files automatically. + +The script can be run with different options to download the different versions of the dataset: + +- `./download_data.py --ds 0` for the full-resolution dataset +- `./download_data.py --ds 1` for the downsampled variant +- `./download_data.py --ds 2` for the two times downsampled variant + +Please move the downloaded datasets to `~/vofflow/data`. + +## Reproducing the Results + +### Prerequisites + +Our reproducibility scripts assume that: +- The ParaView binary is installed at `~/vofflow/install/bin/paraview` (and `pvbatch`, etc.) +- The Plugin is installed at `~/vofflow/install/lib/paraview-5.12/plugins/VofFlow/VofFlow.so` +- The Plugin tools are installed at `~/vofflow/install/bin/SeedGrid` (and `PolyGap`, `SmoothNormals`, `MaxVelocity`) +- The dataset "ds2" was downloaded at `~/vofflow/data/jet-collision-ds2/jet-collision-ds2.pvd` + +All scripts have the paths configured at the top. +If your setup varies, you can update the paths accordingly, i.e., if you use one of the other datasets. + +### Run Separation Boundaries + +This step runs our algorithm and writes all resulting data, i.e., the separation boundaries, to disk in the VTK data formats. +The output is written to `~/vofflow/output/result` + +```shell +~/vofflow/install/bin/pvbatch separation_boundaries.py +``` + +Due to the computation time, we recommend running this step with MPI: +```shell +mpirun -np 32 ~/vofflow/install/bin/pvbatch separation_boundaries.py +``` + +### Run Surface Smoothing Post Processing Step + +```shell +cd ~/vofflow/output/result +~/vofflow/install/bin/SeedGrid seeds.pvd seed_grid.vti +~/vofflow/install/bin/PolyGap bound.pvd seed_grid.vti bound_gap.vtp +~/vofflow/install/bin/SmoothNormals bound_gap.vtp bound_gap_smoothed.vtp +``` + +### Run Separation Boundary Rendering + +Generates one timestep from Figure 14. + +```shell +~/vofflow/install/bin/pvbatch figure_jet_bounds.py +``` + +### Run PLIC Rendering + +Generates the PLIC surfaces of the jet shown in Figure 13. + +```shell +~/vofflow/install/bin/pvbatch figure_jet_plic.py +``` diff --git a/docs/img/plugin_pipeline.png b/docs/img/plugin_pipeline.png new file mode 100644 index 0000000..55f9018 Binary files /dev/null and b/docs/img/plugin_pipeline.png differ diff --git a/docs/img/plugin_properties.png b/docs/img/plugin_properties.png new file mode 100644 index 0000000..f2ee967 Binary files /dev/null and b/docs/img/plugin_properties.png differ diff --git a/docs/img/teaser.png b/docs/img/teaser.png new file mode 100644 index 0000000..c207f76 Binary files /dev/null and b/docs/img/teaser.png differ diff --git a/docs/scripts/build_plugin_linux.sh b/docs/scripts/build_plugin_linux.sh new file mode 100755 index 0000000..1729a90 --- /dev/null +++ b/docs/scripts/build_plugin_linux.sh @@ -0,0 +1,53 @@ +#!/bin/bash + +# This script is based on the ParaViewEasyPluginBuilder +# https://gitlab.kitware.com/paraview/paraview-easy-plugin-builder + +# Create output dir if not exists +if [[ -e plugin_build ]]; then + echo "Error: Directory 'plugin_build' already exists. Please remove!" + exit 1 +fi +mkdir plugin_build + +# Create build script running in the docker container +cat > "plugin_build/run.sh" < lib +mv /plugin_build/lib64 /plugin_build/lib + +# Cleanup file owners +chown -R $(id -u):$(id -g) /plugin_build +EOF +chmod +x "plugin_build/run.sh" + +# Execute build in docker +docker run --rm -v "$(pwd)/plugin_build:/plugin_build" kitware/paraview_org-plugin-devel:5.12.0 /bin/bash -c "/plugin_build/run.sh" + +# Cleanup +rm "plugin_build/run.sh" diff --git a/docs/scripts/download_data.py b/docs/scripts/download_data.py new file mode 100755 index 0000000..9cfed74 --- /dev/null +++ b/docs/scripts/download_data.py @@ -0,0 +1,45 @@ +#!/usr/bin/env python3 + +import argparse +import os +import requests + + +def main(ds: int = 0): + doi = '10.18419/darus-4225' + if ds == 0: + download_dir = 'jet-collision' + elif ds == 1: + download_dir = 'jet-collision-ds1' + elif ds == 2: + download_dir = 'jet-collision-ds2' + else: + raise RuntimeError('Unknown dataset version') + + # Get file list + r = requests.get('https://darus.uni-stuttgart.de/api/datasets/:persistentId/?persistentId=doi:{}'.format(doi)) + dataset_info = r.json() + + # Sort by directory + directories = {} + for f in dataset_info['data']['latestVersion']['files']: + dirname = f['directoryLabel'] + if dirname not in directories: + directories[dirname] = [] + directories[dirname].append((f['dataFile']['id'], f['dataFile']['filename'])) + + # Download + os.makedirs(download_dir, exist_ok=True) + for fileid, filename in directories[download_dir]: + print(filename) + with requests.get('https://darus.uni-stuttgart.de/api/access/datafile/{}'.format(fileid), stream=True) as r, \ + open(os.path.join(download_dir, filename), 'wb') as f: + for chunk in r.iter_content(chunk_size=1024): + f.write(chunk) + + +if __name__ == '__main__': + parser = argparse.ArgumentParser() + parser.add_argument('-d', '--ds', type=int, choices=range(0, 3), default=0) + args = parser.parse_args() + main(args.ds) diff --git a/docs/scripts/figure_jet_bounds.py b/docs/scripts/figure_jet_bounds.py new file mode 100644 index 0000000..cb048a4 --- /dev/null +++ b/docs/scripts/figure_jet_bounds.py @@ -0,0 +1,92 @@ +import os +from paraview.simple import * + +# Config +dataset = '~/vofflow/output/result/bound_gap_smoothed.vtp' +plugin_path = '~/vofflow/install/lib/paraview-5.12/plugins/VofFlow/VofFlow.so' +output_path = '~/vofflow/output/fig_jet_bound.png' + +label_jet = 15 # Use 20 for full dataset, 20 for "ds1", and 15 for "ds2" +label_drop = 8 # Use 12 for full dataset, 12 for "ds1", and 8 for "ds2" + +# Load plugin +LoadPlugin(os.path.expanduser(plugin_path), ns=globals()) + +LoadPalette(paletteName='WhiteBackground') + +renderView1 = GetActiveViewOrCreate('RenderView') + +# Read smoothed bound data +reader = XMLPolyDataReader(registrationName='bounds', FileName=[os.path.expanduser(dataset)]) +reader.PointArrayStatus = ['Labels', 'Normals'] + +# Cutout jet label +threshold1 = Threshold(registrationName='Threshold1', Input=reader) +threshold1.Scalars = ['POINTS', 'Labels'] +threshold1.LowerThreshold = label_jet +threshold1.UpperThreshold = label_jet + +# Mirror for esthetic reasons, that jets flows from left to right in image +reflect1 = Reflect(registrationName='Reflect1', Input=threshold1) +reflect1.Plane = 'Z Min' +reflect1.CopyInput = 0 +reflect1.FlipAllInputArrays = 0 + +# Cut bounds in mirror plane to view "into" them +clip1 = Clip(registrationName='Clip1', Input=reflect1) +clip1.ClipType = 'Plane' +clip1.ClipType.Normal = [0.0, 0.0, 1.0] + +# Cutout drop label +threshold2 = Threshold(registrationName='Threshold2', Input=reader) +threshold2.Scalars = ['POINTS', 'Labels'] +threshold2.LowerThreshold = label_drop +threshold2.UpperThreshold = label_drop + +# Mirror for esthetic reasons, that jets flows from left to right in image +reflect2 = Reflect(registrationName='Reflect2', Input=threshold2) +reflect2.Plane = 'Z Min' +reflect2.CopyInput = 0 +reflect2.FlipAllInputArrays = 0 + +# Cut bounds in mirror plane to view "into" them +clip2 = Clip(registrationName='Clip2', Input=reflect2) +clip2.ClipType = 'Plane' +clip2.ClipType.Normal = [0.0, 0.0, 1.0] + +# Render +color_jet = [1.0, 0.4980392156862745, 0.054901960784313725] +color_drop = [0.12156862745098039, 0.4666666666666667, 0.7058823529411765] + +clip1Display = Show(clip1, renderView1, 'UnstructuredGridRepresentation') +clip1Display.Representation = 'Surface' +clip1Display.ColorArrayName = [None, ''] +clip1Display.AmbientColor = color_jet +clip1Display.DiffuseColor = color_jet +clip1Display.Specular = 0.5 +clip1Display.SpecularPower = 25.0 +clip1Display.SelectNormalArray = 'Normals' + +clip2Display = Show(clip2, renderView1, 'UnstructuredGridRepresentation') +clip2Display.Representation = 'Surface' +clip2Display.ColorArrayName = [None, ''] +clip2Display.AmbientColor = color_drop +clip2Display.DiffuseColor = color_drop +clip2Display.Specular = 0.5 +clip2Display.SpecularPower = 25.0 +clip2Display.SelectNormalArray = 'Normals' + +# Camera setup +renderView1.ViewSize = [2400, 800] +renderView1.ResetActiveCameraToPositiveZ() +renderView1.ResetCamera(False) +renderView1.CameraPosition = [0.29, -0.024, 0.4] +renderView1.CameraFocalPoint = [0.29, -0.024, 0.0] +renderView1.CameraViewUp = [0.0, 1.0, 0.0] +renderView1.OrientationAxesVisibility = 0 +renderView1.Update() + +# Save screenshot +output_path = os.path.expanduser(output_path) +os.makedirs(os.path.dirname(output_path), exist_ok=True) +SaveScreenshot(output_path, renderView1, ImageResolution=[2400, 800]) diff --git a/docs/scripts/figure_jet_plic.py b/docs/scripts/figure_jet_plic.py new file mode 100644 index 0000000..cdb633c --- /dev/null +++ b/docs/scripts/figure_jet_plic.py @@ -0,0 +1,106 @@ +import os +from paraview.simple import * + +# Config +dataset = '~/vofflow/data/jet-collision-ds2/jet-collision-ds2.0152.vtr' +plugin_path = '~/vofflow/install/lib/paraview-5.12/plugins/VofFlow/VofFlow.so' +output_path = '~/vofflow/output/jet_plic.png' + +# Load plugin +LoadPlugin(os.path.expanduser(plugin_path), ns=globals()) + +LoadPalette(paletteName='WhiteBackground') + +renderView1 = GetActiveViewOrCreate('RenderView') + +# Read dataset +reader = XMLRectilinearGridReader(registrationName='dataset', FileName=os.path.expanduser(dataset)) +reader.CellArrayStatus = ['f3-function[-]', 'n_c_3ph[1]', 'vof-function[-]'] + +# PLIC3 Filter +pLIC31 = PLIC3(registrationName='PLIC31', InputGrid=reader) +pLIC31.VoF3 = ['CELLS', 'f3-function[-]'] +pLIC31.VoF = ['CELLS', 'vof-function[-]'] +pLIC31.VoFNorm = ['CELLS', 'n_c_3ph[1]'] +pLIC31.Epsilon = 1e-05 +pLIC31.NumIterations = 15 + +# Filter for surface of inner phase +threshold1 = Threshold(registrationName='Threshold1', Input=pLIC31) +threshold1.Scalars = ['CELLS', 'PhaseType'] +threshold1.LowerThreshold = 1.0 +threshold1.UpperThreshold = 1.0 + +# Filter for surface of outer phase +threshold2 = Threshold(registrationName='Threshold2', Input=pLIC31) +threshold2.Scalars = ['CELLS', 'PhaseType'] +threshold2.LowerThreshold = 2.0 +threshold2.UpperThreshold = 2.0 + +# Mirror the inner surface +reflect1 = Reflect(registrationName='Reflect1', Input=threshold1) +reflect1.Plane = 'Z Min' +reflect1.CopyInput = 0 +reflect1.FlipAllInputArrays = 0 + +# Mirror the outer surface +reflect2 = Reflect(registrationName='Reflect2', Input=threshold2) +reflect2.Plane = 'Z Min' +reflect2.CopyInput = 0 +reflect2.FlipAllInputArrays = 0 + +# Surface colors +color_inner = [0.0, 1.0, 0.0] +color_outer = [0.8784313725490196, 0.8784313725490196, 0.8784313725490196] + +# Display inner phase +threshold1Display = Show(threshold1, renderView1, 'UnstructuredGridRepresentation') +threshold1Display.Representation = 'Surface' +threshold1Display.AmbientColor = color_inner +threshold1Display.ColorArrayName = [None, ''] +threshold1Display.DiffuseColor = color_inner +threshold1Display.Specular = 0.5 +threshold1Display.SpecularPower = 25.0 + +# Display outer phase +threshold2Display = Show(threshold2, renderView1, 'UnstructuredGridRepresentation') +threshold2Display.Representation = 'Surface' +threshold2Display.AmbientColor = color_outer +threshold2Display.ColorArrayName = [None, ''] +threshold2Display.DiffuseColor = color_outer +threshold2Display.Opacity = 0.5 +threshold2Display.Specular = 0.5 +threshold2Display.SpecularPower = 25.0 + +# Display mirrored inner phase +reflect1Display = Show(reflect1, renderView1, 'UnstructuredGridRepresentation') +reflect1Display.Representation = 'Surface' +reflect1Display.AmbientColor = color_inner +reflect1Display.ColorArrayName = [None, ''] +reflect1Display.DiffuseColor = color_inner +reflect1Display.Specular = 0.5 +reflect1Display.SpecularPower = 25.0 + +# Display mirrored outer phase +reflect2Display = Show(reflect2, renderView1, 'UnstructuredGridRepresentation') +reflect2Display.Representation = 'Surface' +reflect2Display.AmbientColor = color_outer +reflect2Display.ColorArrayName = [None, ''] +reflect2Display.DiffuseColor = color_outer +reflect2Display.Opacity = 0.5 +reflect2Display.Specular = 0.5 +reflect2Display.SpecularPower = 25.0 + +# Camera setup +renderView1.ResetActiveCameraToPositiveZ() +renderView1.ResetCamera(False) +renderView1.CameraPosition = [0.34060654044151306, -0.03483813628554344, 0.4401733070104946] +renderView1.CameraFocalPoint = [0.34060654044151306, -0.03483813628554344, 0.005889910257715658] +renderView1.CameraViewUp = [0.0, 1.0, 0.0] +renderView1.OrientationAxesVisibility = 0 +renderView1.Update() + +# Save screenshot +output_path = os.path.expanduser(output_path) +os.makedirs(os.path.dirname(output_path), exist_ok=True) +SaveScreenshot(output_path, renderView1, ImageResolution=[2400, 800]) diff --git a/docs/scripts/requirements.txt b/docs/scripts/requirements.txt new file mode 100644 index 0000000..d80d9fc --- /dev/null +++ b/docs/scripts/requirements.txt @@ -0,0 +1 @@ +requests==2.32.3 diff --git a/docs/scripts/separation_boundaries.py b/docs/scripts/separation_boundaries.py new file mode 100644 index 0000000..19dcb46 --- /dev/null +++ b/docs/scripts/separation_boundaries.py @@ -0,0 +1,63 @@ +import os +from paraview.simple import * + +# Config +dataset = '~/vofflow/data/jet-collision-ds2/jet-collision-ds2.pvd' +plugin_path = '~/vofflow/install/lib/paraview-5.12/plugins/VofFlow/VofFlow.so' +output_path = '~/vofflow/output/result' + +init_timestep = 68 # The time steps used in the paper (fig. 14) are: 68, 82, 96, 110, 124, 138, 152. +ghost_cells = 10 # MPI ghost cells, use 36 for full dataset, 20 for "ds1", and 10 for "ds2". + +# Load plugin +LoadPlugin(os.path.expanduser(plugin_path), ns=globals()) + +# Read dataset +reader = PVDReader(registrationName='dataset', FileName=os.path.expanduser(dataset)) +reader.CellArrays = ['vof-function[-]', 'velocity[cm/s]', 'f3-function[-]', 'n_c_3ph[1]'] +time_steps = reader.TimestepValues + +# Configure VofTracking +vofTracking1 = VofTracking(registrationName='VofTracking1', InputGrid=reader) +vofTracking1.UseThreePhase = 1 +vofTracking1.VoF = ['CELLS', 'vof-function[-]'] +vofTracking1.VoF3 = ['CELLS', 'f3-function[-]'] +vofTracking1.VoFNorm = ['CELLS', 'n_c_3ph[1]'] +vofTracking1.Velocity = ['CELLS', 'velocity[cm/s]'] +vofTracking1.UseComponents = 0 +vofTracking1.ComponentsVoF = ['CELLS', 'None'] +vofTracking1.ComponentsVoF3 = ['CELLS', 'None'] +vofTracking1.UseTargetTimeStep = 1 +vofTracking1.InitTimeStep = init_timestep +vofTracking1.TargetTimeStep = 152 +vofTracking1.Refinement = 2 +vofTracking1.NeighborCorrection = 1 +vofTracking1.CellCorrection = 1 +vofTracking1.PLICCorrection = 1 +vofTracking1.IntegrationMethod = 'RK4' # 'Euler' +vofTracking1.IntegrationSubSteps = 8 +vofTracking1.Epsilon = 1e-05 +vofTracking1.NumIterations = 15 +vofTracking1.GhostCells = ghost_cells +vofTracking1.CutLabels = 0 +vofTracking1.LabelCutType = 'Plane' +vofTracking1.BoundaryMethod = 'DiscreteFlyingEdges3D' # 'DiscreteMarchingCubes' +vofTracking1.OutputDataType = 'vtkImageData' # 'vtkRectilinearGrid' +vofTracking1.OutputState = 1 +vofTracking1.OutputTimeMeasure = 1 +vofTracking1.MirrorXMin = 0 +vofTracking1.MirrorXMax = 0 +vofTracking1.MirrorYMin = 0 +vofTracking1.MirrorYMax = 0 +vofTracking1.MirrorZMin = 1 +vofTracking1.MirrorZMax = 0 + +# Check output dir +output_path = os.path.expanduser(output_path) +os.makedirs(output_path, exist_ok=True) + +# Save output +SaveData(os.path.join(output_path, 'grid.pvd'), proxy=OutputPort(vofTracking1, 0), CompressorType='ZLib', CompressionLevel='1') +SaveData(os.path.join(output_path, 'seeds.pvd'), proxy=OutputPort(vofTracking1, 1), CompressorType='ZLib', CompressionLevel='1') +SaveData(os.path.join(output_path, 'part.pvd'), proxy=OutputPort(vofTracking1, 2), CompressorType='ZLib', CompressionLevel='1') +SaveData(os.path.join(output_path, 'bound.pvd'), proxy=OutputPort(vofTracking1, 3), CompressorType='ZLib', CompressionLevel='1') diff --git a/src/ParaViewPlugin/CMakeLists.txt b/src/ParaViewPlugin/CMakeLists.txt new file mode 100644 index 0000000..0857dfe --- /dev/null +++ b/src/ParaViewPlugin/CMakeLists.txt @@ -0,0 +1,24 @@ +# Overwrite policy version inherited vom ParaView CMake. +cmake_policy(PUSH) +cmake_policy(VERSION 3.12...3.29) + +paraview_add_plugin(VofFlow + VERSION "1.0" + SERVER_MANAGER_XML VofFlow.xml + MODULES + VofFlow::VofDataUtils + VofFlow::VofTracking + VofFlow::Plic + MODULE_FILES + "${CMAKE_CURRENT_SOURCE_DIR}/../VofDataUtils/vtk.module" + "${CMAKE_CURRENT_SOURCE_DIR}/../VofTracking/vtk.module" + "${CMAKE_CURRENT_SOURCE_DIR}/../Plic/vtk.module" + UI_RESOURCES VofFlow.qrc) + +target_link_libraries(VofFlow + PRIVATE + VofFlow::VofDataUtils + VofFlow::VofTracking + VofFlow::Plic) + +cmake_policy(POP) diff --git a/src/ParaViewPlugin/IconPlic.png b/src/ParaViewPlugin/IconPlic.png new file mode 100644 index 0000000..52e1688 Binary files /dev/null and b/src/ParaViewPlugin/IconPlic.png differ diff --git a/src/ParaViewPlugin/IconPlic3.png b/src/ParaViewPlugin/IconPlic3.png new file mode 100644 index 0000000..b1dd18f Binary files /dev/null and b/src/ParaViewPlugin/IconPlic3.png differ diff --git a/src/ParaViewPlugin/IconVofTracking.png b/src/ParaViewPlugin/IconVofTracking.png new file mode 100644 index 0000000..0b1f37d Binary files /dev/null and b/src/ParaViewPlugin/IconVofTracking.png differ diff --git a/src/ParaViewPlugin/VofFlow.qrc b/src/ParaViewPlugin/VofFlow.qrc new file mode 100644 index 0000000..fdfa896 --- /dev/null +++ b/src/ParaViewPlugin/VofFlow.qrc @@ -0,0 +1,7 @@ + + + IconVofTracking.png + IconPlic.png + IconPlic3.png + + diff --git a/src/ParaViewPlugin/VofFlow.xml b/src/ParaViewPlugin/VofFlow.xml new file mode 100644 index 0000000..1ccbe15 --- /dev/null +++ b/src/ParaViewPlugin/VofFlow.xml @@ -0,0 +1,695 @@ + + + + + + + + + + + + + + + + + + + + + + Grid with VoF and velocity arrays. + + + + + + + + + + + + + + + Volume of fluid field. + + + + + + + + + + + + + Volume of fluid field (three-phase). + + + + + + + + + + + + + VoF normals field (three-phase). + + + + + + + + + + Velocity field. + + + + + + + + + + + + + + + + + + Component label field (VoF). + + + + + + + + + + + + + + + + + + Component label field (VoF3). + + + + + + + + + + + + + + + + + + + Number of seed points per cell = (r+1)^3. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + This property specifies the parameters of the cut function + (an implicit description) used to cut the labeled regions. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Grid with VoF array. + + + + + + + + + + + Volume of fluid field. + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Grid with VoF arrays. + + + + + + + + + + + Volume of fluid field (three-phase). + + + + + + + + + + Volume of fluid field. + + + + + + + + + + VoF normals field (three-phase). + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/ParaViewPlugin/paraview.plugin b/src/ParaViewPlugin/paraview.plugin new file mode 100644 index 0000000..14be5ee --- /dev/null +++ b/src/ParaViewPlugin/paraview.plugin @@ -0,0 +1,4 @@ +NAME + VofFlow +DESCRIPTION + VofFlow. diff --git a/src/Plic/CMakeLists.txt b/src/Plic/CMakeLists.txt new file mode 100644 index 0000000..998036c --- /dev/null +++ b/src/Plic/CMakeLists.txt @@ -0,0 +1,13 @@ +# Overwrite policy version inherited vom ParaView CMake. +cmake_policy(PUSH) +cmake_policy(VERSION 3.12...3.29) + +vtk_module_add_module(VofFlow::Plic FORCE_STATIC + CLASSES + vtkAbstractPlic + vtkPlic + vtkPlic3) +target_compile_features(Plic PUBLIC cxx_std_17) +set_target_properties(Plic PROPERTIES CXX_EXTENSIONS OFF) + +cmake_policy(POP) diff --git a/src/Plic/vtk.module b/src/Plic/vtk.module new file mode 100644 index 0000000..4518c1d --- /dev/null +++ b/src/Plic/vtk.module @@ -0,0 +1,6 @@ +NAME + VofFlow::Plic +DEPENDS + VTK::CommonCore + VTK::FiltersCore + VofFlow::VofDataUtils diff --git a/src/Plic/vtkAbstractPlic.cxx b/src/Plic/vtkAbstractPlic.cxx new file mode 100644 index 0000000..9ebcd3b --- /dev/null +++ b/src/Plic/vtkAbstractPlic.cxx @@ -0,0 +1,60 @@ +#include "vtkAbstractPlic.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +vtkStandardNewMacro(vtkAbstractPlic); + +vtkAbstractPlic::vtkAbstractPlic() : Epsilon(0.0), NumIterations(20) { + mpiController_ = vtkMPIController::SafeDownCast(vtkMultiProcessController::GetGlobalController()); +} + +int vtkAbstractPlic::FillInputPortInformation(int vtkNotUsed(port), vtkInformation* info) { + info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkRectilinearGrid"); + return 1; +} + +int vtkAbstractPlic::RequestUpdateExtent(vtkInformation* vtkNotUsed(request), vtkInformationVector** inputVector, + vtkInformationVector* outputVector) { + + vtkInformation* inInfo = inputVector[0]->GetInformationObject(0); + vtkInformation* outInfo = outputVector->GetInformationObject(0); + + int ghostLevels = outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS()); + inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(), std::max(ghostLevels, 1)); + + return 1; +} + +int vtkAbstractPlic::RequestData(vtkInformation* vtkNotUsed(request), vtkInformationVector** inputVector, + vtkInformationVector* outputVector) { + vtkInformation* inInfo = inputVector[0]->GetInformationObject(0); + vtkRectilinearGrid* input = vtkRectilinearGrid::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT())); + if (!input) { + vtkErrorMacro(<< "Input data is missing!"); + return 0; + } + + vtkInformation* outInfo = outputVector->GetInformationObject(0); + vtkPolyData* output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT())); + + try { + return calcPlic(input, output); + } catch (const std::exception& ex) { + vtkErrorMacro(<< ex.what()); + return 0; + } +} + +int vtkAbstractPlic::calcPlic(vtkRectilinearGrid* vtkNotUsed(input), vtkPolyData* vtkNotUsed(output)) { + return 1; +} diff --git a/src/Plic/vtkAbstractPlic.h b/src/Plic/vtkAbstractPlic.h new file mode 100644 index 0000000..2f8b270 --- /dev/null +++ b/src/Plic/vtkAbstractPlic.h @@ -0,0 +1,43 @@ +#pragma once + +#include + +#include "PlicModule.h" + +class vtkMPIController; +class vtkPolyData; +class vtkRectilinearGrid; + +class PLIC_EXPORT vtkAbstractPlic : public vtkPolyDataAlgorithm { +public: + static vtkAbstractPlic* New(); + vtkTypeMacro(vtkAbstractPlic, vtkPolyDataAlgorithm); + + vtkGetMacro(Epsilon, double); + vtkSetMacro(Epsilon, double); + + vtkGetMacro(NumIterations, int); + vtkSetMacro(NumIterations, int); + + vtkAbstractPlic(const vtkAbstractPlic&) = delete; + void operator=(const vtkAbstractPlic&) = delete; + +protected: + vtkAbstractPlic(); + ~vtkAbstractPlic() override = default; + + int FillInputPortInformation(int port, vtkInformation* info) override; + + int RequestUpdateExtent(vtkInformation* request, vtkInformationVector** inputVector, + vtkInformationVector* outputVector) override; + + int RequestData(vtkInformation* request, vtkInformationVector** inputVector, + vtkInformationVector* outputVector) override; + + virtual int calcPlic(vtkRectilinearGrid* input, vtkPolyData* output); + + double Epsilon; + int NumIterations; + + vtkMPIController* mpiController_; +}; diff --git a/src/Plic/vtkPlic.cxx b/src/Plic/vtkPlic.cxx new file mode 100644 index 0000000..54bf562 --- /dev/null +++ b/src/Plic/vtkPlic.cxx @@ -0,0 +1,39 @@ +#include "vtkPlic.h" + +#include +#include +#include + +#include "Grid/DomainInfo.h" +#include "Grid/GridIterator.h" +#include "Grid/GridTypes.h" +#include "Misc/Profiling.h" +#include "Plic/Plic.h" +#include "Plic/PlicPolyData.h" + +vtkStandardNewMacro(vtkPlic); + +int vtkPlic::calcPlic(vtkRectilinearGrid* input, vtkPolyData* output) { + ZoneScoped; + + VofFlow::DomainInfo domainInfo(input, mpiController_); + + vtkDataArray* vofData = GetInputArrayToProcess(0, input); + + VofFlow::vtkPlicPolyData vtkPlic; + + for (const VofFlow::gridCoords_t& e_coords : VofFlow::GridRange(domainInfo.localZeroExtent())) { + const double f = vofData->GetComponent(domainInfo.localExtentCoordToIdx(e_coords), 0); + if (f < Epsilon || f > 1.0 - Epsilon) { + continue; + } + + const auto& plic = VofFlow::calcPlicPoly(domainInfo, domainInfo.localExtentCoordToGridCoord(e_coords), vofData, + Epsilon, NumIterations); + vtkPlic.addPolygon(plic); + } + + output->ShallowCopy(vtkPlic.getPolyData()); + + return 1; +} diff --git a/src/Plic/vtkPlic.h b/src/Plic/vtkPlic.h new file mode 100644 index 0000000..2ed5143 --- /dev/null +++ b/src/Plic/vtkPlic.h @@ -0,0 +1,20 @@ +#pragma once + +#include "vtkAbstractPlic.h" + +#include "PlicModule.h" + +class PLIC_EXPORT vtkPlic : public vtkAbstractPlic { +public: + static vtkPlic* New(); + vtkTypeMacro(vtkPlic, vtkAbstractPlic); + + vtkPlic(const vtkPlic&) = delete; + void operator=(const vtkPlic&) = delete; + +protected: + vtkPlic() = default; + ~vtkPlic() override = default; + + int calcPlic(vtkRectilinearGrid* input, vtkPolyData* output) override; +}; diff --git a/src/Plic/vtkPlic3.cxx b/src/Plic/vtkPlic3.cxx new file mode 100644 index 0000000..a05308d --- /dev/null +++ b/src/Plic/vtkPlic3.cxx @@ -0,0 +1,93 @@ +#include "vtkPlic3.h" + +#include + +#include +#include +#include + +#include "Grid/DomainInfo.h" +#include "Grid/GridIterator.h" +#include "Grid/GridTypes.h" +#include "Misc/CgalUtil.h" +#include "Misc/Profiling.h" +#include "Misc/VofData.h" +#include "Plic/Plic3.h" +#include "Plic/PlicPolyData.h" +#include "Plic/Polygon3D.h" + +vtkStandardNewMacro(vtkPlic3); + +int vtkPlic3::calcPlic(vtkRectilinearGrid* input, vtkPolyData* output) { + ZoneScoped; + + VofFlow::DomainInfo domainInfo(input, mpiController_); + + VofFlow::VofData vofData{ + GetInputArrayToProcess(0, input), + GetInputArrayToProcess(1, input), + GetInputArrayToProcess(2, input), + }; + + VofFlow::vtkPlic3PolyData vtkPlic; + + for (const VofFlow::gridCoords_t& e_coords : VofFlow::GridRange(domainInfo.localZeroExtent())) { + const auto& plic = VofFlow::calcPlic3CellClass(domainInfo, domainInfo.localExtentCoordToGridCoord(e_coords), + vofData, Epsilon, NumIterations); + + if (plic.cellClass == VofFlow::CellClass::EMPTY || plic.cellClass == VofFlow::CellClass::FULL_PHASE1 || + plic.cellClass == VofFlow::CellClass::FULL_PHASE2) { + continue; + } + + // PhaseType: If interface of first phase value is 1, otherwise if interface of second phase, value is 2. + // CellClass must be one of: INTERFACE_PH1_PH2, INTERFACE_PHASE1, INTERFACE_PHASE2, INTERFACE_ALL + const unsigned char phaseType = (plic.cellClass == VofFlow::CellClass::INTERFACE_PHASE2) ? 2 : 1; + + const auto& plic1 = plic.plic1.value(); + vtkPlic.addPolygon(plic.cellClass, phaseType, {plic1.cell.getPolygon(0), plic1.iter, plic1.err}); + + if (plic.cellClass == VofFlow::CellClass::INTERFACE_ALL) { + const auto& plic2 = plic.plic2.value(); + const auto& plane1 = plic1.cell.getPlane(0); + const auto& plane2 = plic2.cell.getPlane(1); + + std::vector polyPoints; + const auto& poly = plic2.cell.getPolygon(1); // Store temporary object + for (const auto& p : poly.points()) { + if (plane1.has_on_positive_side(p)) { + polyPoints.push_back(p); + } + } + + const auto result = CGAL::intersection(plane1, plane2); + if (result) { + if (const VofFlow::K::Line_3* l = boost::get(&*result)) { + const auto result2 = CGAL::intersection(plic1.cell.cellCube(), *l); + if (result2) { + if (const VofFlow::K::Segment_3* s = boost::get(&*result2)) { + polyPoints.push_back(s->source()); + polyPoints.push_back(s->target()); + } else { + const VofFlow::K::Point_3* p = boost::get(&*result2); + polyPoints.push_back(*p); + } + } + } else { + // Planes are identical. No need to add a second plane at the very same position. + } + } else { + // Parallel planes. `has_on_positive_side` has either added all or none points. + } + + if (polyPoints.size() >= 3) { + // Second plane in three-phase cell has always PhaseType 2 + vtkPlic.addPolygon(plic.cellClass, 2, {VofFlow::Polygon3D(plane2, polyPoints), plic2.iter, plic2.err}); + } + } + } + + output->ShallowCopy(vtkPlic.getPolyData()); + + return 1; +} diff --git a/src/Plic/vtkPlic3.h b/src/Plic/vtkPlic3.h new file mode 100644 index 0000000..29e8b5b --- /dev/null +++ b/src/Plic/vtkPlic3.h @@ -0,0 +1,20 @@ +#pragma once + +#include "vtkAbstractPlic.h" + +#include "PlicModule.h" + +class PLIC_EXPORT vtkPlic3 : public vtkAbstractPlic { +public: + static vtkPlic3* New(); + vtkTypeMacro(vtkPlic3, vtkAbstractPlic); + + vtkPlic3(const vtkPlic3&) = delete; + void operator=(const vtkPlic3&) = delete; + +protected: + vtkPlic3() = default; + ~vtkPlic3() override = default; + + int calcPlic(vtkRectilinearGrid* input, vtkPolyData* output) override; +}; diff --git a/src/Tools/CMakeLists.txt b/src/Tools/CMakeLists.txt new file mode 100644 index 0000000..09af4c6 --- /dev/null +++ b/src/Tools/CMakeLists.txt @@ -0,0 +1,66 @@ +# Tools + +find_package(OpenMP REQUIRED) +find_package(nlohmann_json CONFIG REQUIRED) + +add_executable(MaxVelocity MaxVelocity.cpp) +target_compile_features(MaxVelocity PUBLIC cxx_std_17) +set_target_properties(MaxVelocity PROPERTIES CXX_EXTENSIONS OFF) +target_link_libraries(MaxVelocity + PUBLIC + ParaView::VTKExtensionsIOCore + VofDataUtils + OpenMP::OpenMP_CXX) +if (UNIX) + set_target_properties(MaxVelocity PROPERTIES + BUILD_RPATH_USE_ORIGIN TRUE + INSTALL_RPATH "\$ORIGIN/../lib") +endif () + +add_executable(SeedGrid SeedGrid.cpp) +target_compile_features(SeedGrid PUBLIC cxx_std_17) +set_target_properties(SeedGrid PROPERTIES CXX_EXTENSIONS OFF) +target_link_libraries(SeedGrid + PUBLIC + ParaView::VTKExtensionsIOCore + nlohmann_json::nlohmann_json + VofDataUtils + VofTracking) +if (UNIX) + set_target_properties(SeedGrid PROPERTIES + BUILD_RPATH_USE_ORIGIN TRUE + INSTALL_RPATH "\$ORIGIN/../lib") +endif () + +add_executable(PolyGap PolyGap.cpp) +target_compile_features(PolyGap PUBLIC cxx_std_17) +set_target_properties(PolyGap PROPERTIES CXX_EXTENSIONS OFF) +target_link_libraries(PolyGap + PUBLIC + ParaView::VTKExtensionsIOCore + VofDataUtils) +if (UNIX) + set_target_properties(PolyGap PROPERTIES + BUILD_RPATH_USE_ORIGIN TRUE + INSTALL_RPATH "\$ORIGIN/../lib") +endif () + +add_executable(SmoothNormals SmoothNormals.cpp) +target_compile_features(SmoothNormals PUBLIC cxx_std_17) +set_target_properties(SmoothNormals PROPERTIES CXX_EXTENSIONS OFF) +target_link_libraries(SmoothNormals + PUBLIC + ParaView::VTKExtensionsIOCore + VTK::FiltersParallel + VofDataUtils) +if (UNIX) + set_target_properties(SmoothNormals PROPERTIES + BUILD_RPATH_USE_ORIGIN TRUE + INSTALL_RPATH "\$ORIGIN/../lib") +endif () + +# Install +include(GNUInstallDirs) + +install(TARGETS MaxVelocity SeedGrid PolyGap SmoothNormals + RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}) diff --git a/src/Tools/MaxVelocity.cpp b/src/Tools/MaxVelocity.cpp new file mode 100644 index 0000000..b3cac65 --- /dev/null +++ b/src/Tools/MaxVelocity.cpp @@ -0,0 +1,220 @@ +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Grid/GridTypes.h" +#include "Math/Vector.h" + +namespace fs = std::filesystem; + +static const std::string velName = "velocity[cm/s]"; +static const std::string vof1Name = "vof-function[-]"; +static const std::string vof2Name = "f3-function[-]"; + +struct Result { + explicit Result(float init) : maxAll(init), maxVof1(init), maxVof2(init), maxEmpty(init) {} + + VofFlow::vec3 maxAll; + VofFlow::vec3 maxVof1; + VofFlow::vec3 maxVof2; + VofFlow::vec3 maxEmpty; + + void maxMerge(const Result& other) { + maxAll = VofFlow::max(maxAll, other.maxAll); + maxVof1 = VofFlow::max(maxVof1, other.maxVof1); + maxVof2 = VofFlow::max(maxVof2, other.maxVof2); + maxEmpty = VofFlow::max(maxEmpty, other.maxEmpty); + } +}; + +struct MaxVelocityWorker { + float eps_ = 10e-6; + + Result result; + + MaxVelocityWorker() : result(std::numeric_limits::lowest()) {} + + template + void operator()(vtkAOSDataArrayTemplate* velArray, vtkAOSDataArrayTemplate* vof1Array, + vtkAOSDataArrayTemplate* vof2Array) { + VTK_ASSUME(velArray->GetNumberOfComponents() == 3); + VTK_ASSUME(vof1Array->GetNumberOfComponents() == 1); + VTK_ASSUME(vof2Array->GetNumberOfComponents() == 1); + + if (velArray->GetNumberOfComponents() != 3 || vof1Array->GetNumberOfComponents() != 1 || + vof2Array->GetNumberOfComponents() != 1) { + throw std::runtime_error("Bad array components!"); + } + + vtkIdType num = velArray->GetNumberOfTuples(); + if (num != vof1Array->GetNumberOfTuples() || num != vof2Array->GetNumberOfTuples()) { + throw std::runtime_error("Bad array size!"); + } + + vtkDataArrayAccessor> vel(velArray); + vtkDataArrayAccessor> vof1(vof1Array); + vtkDataArrayAccessor> vof2(vof2Array); + +// MSVC only supports OpenMP 2.0 +#ifndef _WIN32 +#pragma omp declare reduction(max:Result : omp_out.maxMerge(omp_in)) \ + initializer(omp_priv = Result(std::numeric_limits::lowest())) + +#pragma omp parallel for reduction(max : result) +#endif + for (vtkIdType i = 0; i < num; i++) { + const bool isVof1 = vof1.Get(i, 0) >= eps_; + const bool isVof2 = vof2.Get(i, 0) >= eps_; + const bool isEmpty = !(isVof1 || isVof2); + + VofFlow::vec3 v{ + std::abs(static_cast(vel.Get(i, 0))), + std::abs(static_cast(vel.Get(i, 1))), + std::abs(static_cast(vel.Get(i, 2))), + }; + + result.maxAll = VofFlow::max(result.maxAll, v); + if (isVof1) { + result.maxVof1 = VofFlow::max(result.maxVof1, v); + } + if (isVof2) { + result.maxVof2 = VofFlow::max(result.maxVof2, v); + } + if (isEmpty) { + result.maxEmpty = VofFlow::max(result.maxEmpty, v); + } + } + } +}; + +void maxVelocity(vtkSmartPointer& grid, double dt, Result& r) { + auto cellData = grid->GetCellData(); + if (!cellData->HasArray(velName.c_str()) || !cellData->HasArray(vof1Name.c_str()) || + !cellData->HasArray(vof2Name.c_str())) { + throw std::runtime_error("Missing velocity array!"); + } + vtkSmartPointer velArray = cellData->GetArray(velName.c_str()); + vtkSmartPointer vof1Array = cellData->GetArray(vof1Name.c_str()); + vtkSmartPointer vof2Array = cellData->GetArray(vof2Name.c_str()); + + MaxVelocityWorker worker; + + typedef vtkArrayDispatch::Dispatch3ByValueType + Dispatcher; + if (!Dispatcher::Execute(velArray, vof1Array, vof2Array, worker)) { + throw std::runtime_error("Cannot dispatch array worker!"); + } + + VofFlow::extent_t extent; + VofFlow::bounds_t bounds; + grid->GetExtent(extent.data()); + grid->GetBounds(bounds.data()); + auto dims = VofFlow::Grid::extentDimensions(extent); + VofFlow::vec3 cellDimsUniform{ + static_cast((bounds[1] - bounds[0]) / dims[0]), + static_cast((bounds[3] - bounds[2]) / dims[1]), + static_cast((bounds[5] - bounds[4]) / dims[2]), + }; + + // Calc result in cell lengths (assuming uniform). + worker.result.maxAll = worker.result.maxAll * static_cast(dt) / cellDimsUniform; + worker.result.maxVof1 = worker.result.maxVof1 * static_cast(dt) / cellDimsUniform; + worker.result.maxVof2 = worker.result.maxVof2 * static_cast(dt) / cellDimsUniform; + worker.result.maxEmpty = worker.result.maxEmpty * static_cast(dt) / cellDimsUniform; + + r.maxMerge(worker.result); +} + +int main(int argc, char* argv[]) { + if (argc != 2) { + std::cout << "Usage: MaxVelocity " << std::endl; + return 1; + } + fs::path path_data(argv[1]); + if (!fs::is_regular_file(path_data)) { + std::cout << "Input file " << path_data.string() << " does not exist!" << std::endl; + return 1; + } + + // Read data + vtkSmartPointer reader = vtkSmartPointer::New(); + reader->SetFileName(path_data.string().c_str()); + // DisableAllArrays() seems to be required twice to not load data on first Update() + reader->UpdateInformation(); + reader->GetPointDataArraySelection()->DisableAllArrays(); + reader->GetCellDataArraySelection()->DisableAllArrays(); + reader->GetColumnArraySelection()->DisableAllArrays(); + reader->Update(); + reader->GetPointDataArraySelection()->DisableAllArrays(); + reader->GetCellDataArraySelection()->DisableAllArrays(); + reader->GetColumnArraySelection()->DisableAllArrays(); + + // Get number of time steps + int timeFrom, timeTo = 0; + reader->GetTimeStepRange(timeFrom, timeTo); + + // Read time data + std::vector timesteps; + for (int t = timeFrom; t <= timeTo; t++) { + reader->SetTimeStep(t); + reader->Update(); + auto timeArray = reader->GetOutputAsDataSet()->GetFieldData()->GetArray("TimeValue"); + timesteps.push_back(timeArray->GetTuple1(0)); + } + + // Enable required arrays + reader->SetTimeStep(0); + auto cellArr = reader->GetCellDataArraySelection(); + if (!cellArr->ArrayExists(velName.c_str()) || !cellArr->ArrayExists(vof1Name.c_str()) || + !cellArr->ArrayExists(vof2Name.c_str())) { + throw std::runtime_error("Missing arrays!"); + } + cellArr->EnableArray(velName.c_str()); + cellArr->EnableArray(vof1Name.c_str()); + cellArr->EnableArray(vof2Name.c_str()); + + // Calc max velocity + Result r(std::numeric_limits::lowest()); + + for (int t = timeFrom; t <= timeTo; t++) { + std::cout << "Load timestep " << t << "/" << timeTo << std::endl; + const int last_t = std::max(timeFrom, t - 1); + const int next_t = std::min(timeTo, t + 1); + const double dt = std::max(timesteps[t] - timesteps[last_t], timesteps[next_t] - timesteps[t]); + + reader->SetTimeStep(t); + reader->Update(); + vtkSmartPointer grid = vtkRectilinearGrid::SafeDownCast(reader->GetOutputAsDataSet()); + if (grid == nullptr) { + throw std::runtime_error("Cannot load grid!"); + } + + maxVelocity(grid, dt, r); + } + + std::cout << "Max velocity in number of cells (assuming uniform grid) by phase:" << std::endl; + std::cout << "All: " << r.maxAll << std::endl; + std::cout << "Vof1: " << r.maxVof1 << std::endl; + std::cout << "Vof2: " << r.maxVof2 << std::endl; + std::cout << "Empty: " << r.maxEmpty << std::endl; + + return 0; +} diff --git a/src/Tools/PolyGap.cpp b/src/Tools/PolyGap.cpp new file mode 100644 index 0000000..ed77b81 --- /dev/null +++ b/src/Tools/PolyGap.cpp @@ -0,0 +1,163 @@ +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Grid/GridTypes.h" +#include "VtkUtil/VtkUtilWriter.h" + +namespace fs = std::filesystem; + +int main(int argc, char* argv[]) { + if (argc != 4) { + std::cout << "Usage: PolyGap " << std::endl; + return 1; + } + fs::path path_poly(argv[1]); + fs::path path_grid(argv[2]); + fs::path path_out(argv[3]); + if (!fs::is_regular_file(path_poly)) { + std::cout << "Input file " << path_poly.string() << " does not exist!" << std::endl; + return 1; + } + if (!fs::is_regular_file(path_grid)) { + std::cout << "Input file " << path_grid.string() << " does not exist!" << std::endl; + return 1; + } + if (fs::exists(path_out)) { + std::cout << "Output file " << path_out.string() << " already existing!" << std::endl; + return 1; + } + + // Read data + vtkSmartPointer reader = vtkSmartPointer::New(); + reader->SetFileName(path_poly.string().c_str()); + reader->Update(); + vtkSmartPointer bound = vtkPolyData::SafeDownCast(reader->GetOutputAsDataSet()); + + vtkSmartPointer imgReader = vtkSmartPointer::New(); + imgReader->SetFileName(path_grid.string().c_str()); + imgReader->Update(); + vtkSmartPointer grid = imgReader->GetOutput(); + + auto grid_labels = vtkIntArray::SafeDownCast(grid->GetPointData()->GetArray("Labels")); + auto bound_labels = vtkIntArray::SafeDownCast(bound->GetPointData()->GetArray("Labels")); + + VofFlow::posCoords_t spacing; + grid->GetSpacing(spacing.data()); + std::cout << "Spacing: " << spacing[0] << " " << spacing[1] << " " << spacing[2] << std::endl; + + const double displacement = 0.02 * spacing[0]; // TODO assumes uniform for all dimensions + constexpr double eps = 1e-3; + constexpr auto roundCheck = [eps](double& v) { + if (v >= 0.0 && v < eps) { + v = 0.0; + } else if (v > 0.5 - eps && v < 0.5 + eps) { + v = 0.5; + } else if (v > 1.0 - eps && v <= 1.0) { + v = 1.0; + } else { + throw std::runtime_error("Bad value: " + std::to_string(v)); + } + }; + + for (vtkIdType i = 0; i < bound->GetNumberOfPoints(); i++) { + VofFlow::posCoords_t p{}; + bound->GetPoints()->GetPoint(i, p.data()); + + VofFlow::gridCoords_t ijk{}; + VofFlow::posCoords_t pcoords{}; + int result = grid->ComputeStructuredCoordinates(p.data(), ijk.data(), pcoords.data()); + if (result != 1) { + throw std::runtime_error("Out of bounds!"); + } + roundCheck(pcoords[0]); + roundCheck(pcoords[1]); + roundCheck(pcoords[2]); + + // Each point should be in the middle of one line connecting cell edges. + // Therefore, they should have exactly one 0.5 pcoord. + int count_half = 0; + if (pcoords[0] == 0.5) { + count_half++; + } + if (pcoords[1] == 0.5) { + count_half++; + } + if (pcoords[2] == 0.5) { + count_half++; + } + if (count_half != 1) { + throw std::runtime_error("Invalid point setup!"); + } + + // pcoord = 1 is the same as pcoord = 0 with increasing ikj by 1. + // Warning maybe dangerous at upper bound without check, but generated grid has enough margin. + if (pcoords[0] == 1.0) { + pcoords[0] = 0.0; + ijk[0]++; + } + if (pcoords[1] == 1.0) { + pcoords[1] = 0.0; + ijk[1]++; + } + if (pcoords[2] == 1.0) { + pcoords[2] = 0.0; + ijk[2]++; + } + + // Do displacement + const auto b_label = bound_labels->GetValue(i); + if (pcoords[0] == 0.5) { + auto label_lower = grid_labels->GetValue(grid->ComputePointId(ijk.data())); + ijk[0]++; + auto label_upper = grid_labels->GetValue(grid->ComputePointId(ijk.data())); + if (label_lower == b_label && label_upper != b_label) { + p[0] -= displacement; + } else if (label_upper == b_label && label_lower != b_label) { + p[0] += displacement; + } else { + throw std::runtime_error("Invalid labeling found!"); + } + } + if (pcoords[1] == 0.5) { + auto label_lower = grid_labels->GetValue(grid->ComputePointId(ijk.data())); + ijk[1]++; + auto label_upper = grid_labels->GetValue(grid->ComputePointId(ijk.data())); + if (label_lower == b_label && label_upper != b_label) { + p[1] -= displacement; + } else if (label_upper == b_label && label_lower != b_label) { + p[1] += displacement; + } else { + throw std::runtime_error("Invalid labeling found!"); + } + } + if (pcoords[2] == 0.5) { + auto label_lower = grid_labels->GetValue(grid->ComputePointId(ijk.data())); + ijk[2]++; + auto label_upper = grid_labels->GetValue(grid->ComputePointId(ijk.data())); + if (label_lower == b_label && label_upper != b_label) { + p[2] -= displacement; + } else if (label_upper == b_label && label_lower != b_label) { + p[2] += displacement; + } else { + throw std::runtime_error("Invalid labeling found!"); + } + } + + bound->GetPoints()->SetPoint(i, p.data()); + } + + VofFlow::writeData(bound, path_out.string(), true); + + return 0; +} diff --git a/src/Tools/SeedGrid.cpp b/src/Tools/SeedGrid.cpp new file mode 100644 index 0000000..d39121d --- /dev/null +++ b/src/Tools/SeedGrid.cpp @@ -0,0 +1,143 @@ +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Constants.h" +#include "Grid/GridTypes.h" +#include "VtkUtil/VtkUtilWriter.h" + +namespace fs = std::filesystem; + +VofFlow::gridCoords_t idxToSeedCoord(VofFlow::dim_t numPoints, vtkIdType idx) { + const int x = static_cast(idx % numPoints[0]); + const int y = static_cast((idx / numPoints[0]) % numPoints[1]); + const int z = static_cast(idx / (numPoints[0] * numPoints[1])); + return {x, y, z}; +} + +int main(int argc, char* argv[]) { + if (argc != 3) { + std::cout << "Usage: SeedGrid " << std::endl; + return 1; + } + fs::path path_in(argv[1]); + fs::path path_out(argv[2]); + if (!fs::is_regular_file(path_in)) { + std::cout << "Input file " << path_in.string() << " does not exist!" << std::endl; + return 1; + } + if (fs::exists(path_out)) { + std::cout << "Output file " << path_out.string() << " already existing!" << std::endl; + return 1; + } + + // Read dataset + vtkSmartPointer reader = vtkSmartPointer::New(); + reader->SetFileName(path_in.string().c_str()); + reader->Update(); + vtkSmartPointer dataset = vtkPolyData::SafeDownCast(reader->GetOutputAsDataSet()); + if (dataset == nullptr) { + std::cout << "Cannot read poly data!" << std::endl; + return 1; + } + + // Read state + auto state_str = vtkStringArray::SafeDownCast(dataset->GetFieldData()->GetAbstractArray("State"))->GetValue(0); + auto state = nlohmann::json::parse(state_str); + + auto labels = vtkIntArray::SafeDownCast(dataset->GetPointData()->GetArray("Labels")); + auto seed_idx = vtkIdTypeArray::SafeDownCast(dataset->GetPointData()->GetArray("SeedIdx")); + + auto ext = state["domain"]["GlobalExtent"].get(); + auto bounds = state["domain"]["GlobalBounds"].get(); + auto numPointsPerDim = std::max(0, state["parameters"]["Refinement"].get()) + 1; + + VofFlow::dim_t numPoints{ + (ext[1] - ext[0]) * numPointsPerDim, + (ext[3] - ext[2]) * numPointsPerDim, + (ext[5] - ext[4]) * numPointsPerDim, + }; + + VofFlow::dim_t image_dim{numPoints[0] + 4, numPoints[1] + 4, numPoints[2] + 4}; + + VofFlow::dim_t image_min{ + ext[0] * numPointsPerDim - 2, + ext[2] * numPointsPerDim - 2, + ext[4] * numPointsPerDim - 2, + }; + + VofFlow::posCoords_t spacing{ + (bounds[1] - bounds[0]) / numPoints[0], + (bounds[3] - bounds[2]) / numPoints[1], + (bounds[5] - bounds[4]) / numPoints[2], + }; + + VofFlow::posCoords_t origin{ + bounds[0] + (image_min[0] + 0.5) * spacing[0], + bounds[2] + (image_min[1] + 0.5) * spacing[1], + bounds[4] + (image_min[2] + 0.5) * spacing[2], + }; + + // Output image + vtkSmartPointer imgdata = vtkSmartPointer::New(); + imgdata->SetDimensions(image_dim.data()); + imgdata->SetOrigin(origin.data()); + imgdata->SetSpacing(spacing.data()); + imgdata->AllocateScalars(VTK_INT, 1); + imgdata->GetPointData()->GetScalars()->SetName("Labels"); + imgdata->GetPointData()->GetScalars()->Fill(VofFlow::ErrorLabels::MAX_ERROR_LABEL); + + // Add labels + int maxLabel = -1; + + for (int i = 0; i < seed_idx->GetNumberOfTuples(); i++) { + const auto s = idxToSeedCoord(numPoints, seed_idx->GetValue(i)); + const auto idx = imgdata->GetScalarIndex(s[0] - image_min[0], s[1] - image_min[1], s[2] - image_min[2]); + int label = labels->GetValue(i); + imgdata->GetPointData()->GetScalars()->SetTuple1(idx, label); + if (label > maxLabel) { + maxLabel = label; + } + } + std::cout << "Max label: " << maxLabel << std::endl; + + // Calc crop + int min_sx = INT_MAX; + int max_sx = INT_MIN; + int min_sy = INT_MAX; + int max_sy = INT_MIN; + int min_sz = INT_MAX; + int max_sz = INT_MIN; + for (int i = 0; i < seed_idx->GetNumberOfTuples(); i++) { + const auto s = idxToSeedCoord(numPoints, seed_idx->GetValue(i)); + min_sx = std::min(s[0], min_sx); + min_sy = std::min(s[1], min_sy); + min_sz = std::min(s[2], min_sz); + max_sx = std::max(s[0], max_sx); + max_sy = std::max(s[1], max_sy); + max_sz = std::max(s[2], max_sz); + } + + std::cout << "BBox Min: [" << min_sx << " " << min_sy << " " << min_sz << "] Max: [" << max_sx << " " << max_sy + << " " << max_sz << "]" << std::endl; + + VofFlow::extent_t crop{min_sx - 2, max_sx + 4 + 2, min_sy - 2, max_sy + 4 + 2, min_sz - 2, max_sz + 4 + 2}; + imgdata->Crop(crop.data()); + + VofFlow::writeData(imgdata, path_out.string(), true); + + return 0; +} diff --git a/src/Tools/SmoothNormals.cpp b/src/Tools/SmoothNormals.cpp new file mode 100644 index 0000000..2749e5a --- /dev/null +++ b/src/Tools/SmoothNormals.cpp @@ -0,0 +1,67 @@ +#include +#include + +#include +#include +#include +#include +#include + +#include "VtkUtil/VtkUtilWriter.h" + +namespace fs = std::filesystem; + +int main(int argc, char* argv[]) { + if (argc != 3) { + std::cout << "Usage: SmoothNormals " << std::endl; + return 1; + } + fs::path path_in(argv[1]); + fs::path path_out(argv[2]); + if (!fs::is_regular_file(path_in)) { + std::cout << "Input file " << path_in.string() << " does not exist!" << std::endl; + return 1; + } + if (fs::exists(path_out)) { + std::cout << "Output file " << path_out.string() << " already existing!" << std::endl; + return 1; + } + + // Read dataset + vtkSmartPointer reader = vtkSmartPointer::New(); + reader->SetFileName(path_in.string().c_str()); + reader->Update(); + vtkSmartPointer dataset = reader->GetOutput(); + if (dataset == nullptr) { + std::cout << "Cannot read poly data!" << std::endl; + return 1; + } + + // Smooth + vtkSmartPointer smooth = vtkSmartPointer::New(); + smooth->SetNumberOfIterations(200); + smooth->SetConvergence(0.0); + smooth->SetInputData(dataset); + smooth->Update(); + + vtkSmartPointer surfnorm = vtkSmartPointer::New(); + surfnorm->SetFeatureAngle(30.0); + surfnorm->SetSplitting(true); + surfnorm->SetConsistency(true); + surfnorm->SetFlipNormals(false); + surfnorm->SetNonManifoldTraversal(true); + surfnorm->SetComputeCellNormals(false); + surfnorm->SetPieceInvariant(true); + surfnorm->SetInputData(smooth->GetOutput()); + surfnorm->Update(); + + // Replace normals + auto surf_normals = surfnorm->GetOutput()->GetPointData()->GetArray("Normals"); + dataset->GetPointData()->RemoveArray("Normals"); + int idx = dataset->GetPointData()->AddArray(surf_normals); + dataset->GetPointData()->SetActiveAttribute(idx, vtkDataSetAttributes::NORMALS); + + VofFlow::writeData(dataset, path_out.string(), true); + + return 0; +} diff --git a/src/VofDataUtils/CMakeLists.txt b/src/VofDataUtils/CMakeLists.txt new file mode 100644 index 0000000..9f47d9c --- /dev/null +++ b/src/VofDataUtils/CMakeLists.txt @@ -0,0 +1,71 @@ +# Overwrite policy version inherited vom ParaView CMake. +cmake_policy(PUSH) +cmake_policy(VERSION 3.12...3.29) + +find_package(CGAL REQUIRED) + +vtk_module_add_module(VofFlow::VofDataUtils FORCE_STATIC + NOWRAP_HEADERS + Grid/DataInterpolation.h + Grid/DomainInfo.h + Grid/Gradient.h + Grid/GridData.h + Grid/GridIterator.h + Grid/GridTypes.h + Math/Vector.h + Misc/CgalUtil.h + Misc/ListSearch.h + Misc/Profiling.h + Misc/VofData.h + Plic/CachedPlic.h + Plic/Plic.h + Plic/Plic3.h + Plic/PlicPolyData.h + Plic/PlicUtil.h + Plic/PolyCell.h + Plic/Polygon3D.h + VtkUtil/VtkUtilArray.h + VtkUtil/VtkUtilMPI.h + VtkUtil/VtkUtilWriter.h + SOURCES + Grid/DataInterpolation.cpp + Grid/DomainInfo.cpp + Grid/Gradient.cpp + Plic/Plic.cpp + Plic/Plic3.cpp + Plic/PlicPolyData.cpp + Plic/PlicUtil.cpp + Plic/PolyCell.cpp + Plic/Polygon3D.cpp + VtkUtil/VtkUtilArray.cpp + VtkUtil/VtkUtilWriter.cpp) +target_compile_features(VofDataUtils PUBLIC cxx_std_17) +set_target_properties(VofDataUtils PROPERTIES CXX_EXTENSIONS OFF) + +target_link_libraries(VofDataUtils PUBLIC CGAL::CGAL) + +option(VOFFLOW_USE_TRACY "" OFF) +if (VOFFLOW_USE_TRACY) + # We require TRACY_DELAYED_INIT. This would require an vcpkg overlay port, use fetch content instead. + FetchContent_Declare(tracy + URL "https://github.com/wolfpld/tracy/archive/v0.10.tar.gz" + URL_HASH SHA256=a76017d928f3f2727540fb950edd3b736caa97b12dbb4e5edce66542cbea6600) + FetchContent_GetProperties(tracy) + if (NOT tracy_POPULATED) + message(STATUS "Fetch tracy ...") + FetchContent_Populate(tracy) + option(TRACY_ENABLE "" ON) + option(TRACY_DELAYED_INIT "" OFF) + option(TRACY_STATIC "" OFF) + add_subdirectory(${tracy_SOURCE_DIR} ${tracy_BINARY_DIR}) # TODO EXCLUDE_FROM_ALL but need to install .so + mark_as_advanced(FORCE + FETCHCONTENT_SOURCE_DIR_TRACY + FETCHCONTENT_UPDATES_DISCONNECTED_TRACY) + endif () + + target_link_libraries(VofDataUtils PUBLIC TracyClient) + + target_compile_definitions(VofDataUtils PUBLIC VOFFLOW_USE_TRACY) +endif () + +cmake_policy(POP) diff --git a/src/VofDataUtils/Grid/DataInterpolation.cpp b/src/VofDataUtils/Grid/DataInterpolation.cpp new file mode 100644 index 0000000..69eb6e6 --- /dev/null +++ b/src/VofDataUtils/Grid/DataInterpolation.cpp @@ -0,0 +1,108 @@ +#include "DataInterpolation.h" + +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace { + struct CellDataInterpolation3Worker { + + const VofFlow::DomainInfo& domainInfo_; + const VofFlow::gridCoords_t& cellCoords_; + const VofFlow::posCoords_t& pCoords_; + VofFlow::vec3 result_; + + CellDataInterpolation3Worker(const VofFlow::DomainInfo& domainInfo, const VofFlow::gridCoords_t& cellCoords, + const VofFlow::posCoords_t& pCoords) + : domainInfo_(domainInfo), + cellCoords_(cellCoords), + pCoords_(pCoords), + result_(0.0f, 0.0f, 0.0f) {} + + template + void operator()(vtkAOSDataArrayTemplate* dataArray) { + VTK_ASSUME(dataArray->GetNumberOfComponents() == 3); + + vtkDataArrayAccessor> data(dataArray); + + int lx = cellCoords_[0]; + int ly = cellCoords_[1]; + int lz = cellCoords_[2]; + + double x = pCoords_[0] - 0.5; + double y = pCoords_[1] - 0.5; + double z = pCoords_[2] - 0.5; + + if (pCoords_[0] < 0.5) { + lx -= 1; + x = pCoords_[0] + 0.5; + } + if (pCoords_[1] < 0.5) { + ly -= 1; + y = pCoords_[1] + 0.5; + } + if (pCoords_[2] < 0.5) { + lz -= 1; + z = pCoords_[2] + 0.5; + } + + int ux = lx + 1; + int uy = ly + 1; + int uz = lz + 1; + + lx = std::max(lx, 0); + ly = std::max(ly, 0); + lz = std::max(lz, 0); + ux = std::min(ux, domainInfo_.cellDims()[0] - 1); + uy = std::min(uy, domainInfo_.cellDims()[1] - 1); + uz = std::min(uz, domainInfo_.cellDims()[2] - 1); + + const std::array idx{ + domainInfo_.gridCoordToIdx(lx, ly, lz), + domainInfo_.gridCoordToIdx(ux, ly, lz), + domainInfo_.gridCoordToIdx(lx, uy, lz), + domainInfo_.gridCoordToIdx(ux, uy, lz), + domainInfo_.gridCoordToIdx(lx, ly, uz), + domainInfo_.gridCoordToIdx(ux, ly, uz), + domainInfo_.gridCoordToIdx(lx, uy, uz), + domainInfo_.gridCoordToIdx(ux, uy, uz), + }; + std::array values; + for (int i = 0; i < 8; i++) { + values[i] = VofFlow::vec3{ + static_cast(data.Get(idx[i], 0)), + static_cast(data.Get(idx[i], 1)), + static_cast(data.Get(idx[i], 2)), + }; + } + + const VofFlow::vec3 a = VofFlow::lerp(values[0], values[1], static_cast(x)); + const VofFlow::vec3 b = VofFlow::lerp(values[2], values[3], static_cast(x)); + const VofFlow::vec3 c = VofFlow::lerp(values[4], values[5], static_cast(x)); + const VofFlow::vec3 d = VofFlow::lerp(values[6], values[7], static_cast(x)); + + const VofFlow::vec3 e = VofFlow::lerp(a, b, static_cast(y)); + const VofFlow::vec3 f = VofFlow::lerp(c, d, static_cast(y)); + + result_ = VofFlow::lerp(e, f, static_cast(z)); + } + }; +} // namespace + +VofFlow::vec3 VofFlow::interpolateCellDataVec3(const vtkSmartPointer& data, const DomainInfo& domainInfo, + const gridCoords_t& cellCoords, const posCoords_t& pCoords) { + CellDataInterpolation3Worker worker(domainInfo, cellCoords, pCoords); + + typedef vtkArrayDispatch::DispatchByValueType Dispatcher; + if (!Dispatcher::Execute(data, worker)) { + throw std::runtime_error("Cannot dispatch array worker!"); + } + + return worker.result_; +} diff --git a/src/VofDataUtils/Grid/DataInterpolation.h b/src/VofDataUtils/Grid/DataInterpolation.h new file mode 100644 index 0000000..d8d65ea --- /dev/null +++ b/src/VofDataUtils/Grid/DataInterpolation.h @@ -0,0 +1,26 @@ +#pragma once + +#include +#include + +#include "../Math/Vector.h" +#include "DomainInfo.h" +#include "GridTypes.h" + +namespace VofFlow { + vec3 interpolateCellDataVec3(const vtkSmartPointer& data, const DomainInfo& domainInfo, + const gridCoords_t& cellCoords, const posCoords_t& pCoords); + + inline vec3 interpolateCellDataVec3(const vtkSmartPointer& data, const DomainInfo& domainInfo, + const StructuredCoordinates& coords) { + return interpolateCellDataVec3(data, domainInfo, coords.cellCoords, coords.relativeCoords); + } + + inline vec3 interpolateMultiCellDataVec3(const vtkSmartPointer& data0, + const vtkSmartPointer& data1, const DomainInfo& domainInfo, const vec3& pos, float t) { + const auto coords = domainInfo.computeNearestStructuredCoordinates(pos); + const vec3 v0 = interpolateCellDataVec3(data0, domainInfo, coords); + const vec3 v1 = interpolateCellDataVec3(data1, domainInfo, coords); + return lerp(v0, v1, t); + } +} // namespace VofFlow diff --git a/src/VofDataUtils/Grid/DomainInfo.cpp b/src/VofDataUtils/Grid/DomainInfo.cpp new file mode 100644 index 0000000..79458f1 --- /dev/null +++ b/src/VofDataUtils/Grid/DomainInfo.cpp @@ -0,0 +1,248 @@ +#include "DomainInfo.h" + +#include + +#include +#include +#include +#include +#include +#include + +#include "VtkUtil/VtkUtilArray.h" +#include "VtkUtil/VtkUtilMPI.h" + +namespace { + template + std::vector calcCellSizes(const std::vector& data) { + if (data.size() < 1) { + throw std::runtime_error("Bad array size!"); + } + std::vector result(data.size() - 1); + for (std::size_t i = 0; i < data.size() - 1; i++) { + result[i] = data[i + 1] - data[i]; + } + return result; + } + + template + std::vector calcCellCenters(const std::vector& data) { + if (data.size() < 1) { + throw std::runtime_error("Bad array size!"); + } + std::vector result(data.size() - 1); + for (std::size_t i = 0; i < data.size() - 1; i++) { + result[i] = static_cast(0.5) * (data[i] + data[i + 1]); + } + return result; + } + + template + bool isUniformCoords(const std::vector& data, T eps = std::numeric_limits::epsilon()) { + if (data.size() < 2) { + return true; + } + T min_dist = std::numeric_limits::max(); + T max_dist = std::numeric_limits::lowest(); + for (std::size_t i = 1; i < data.size(); i++) { + const T dist = data[i] - data[i - 1]; + min_dist = std::min(min_dist, dist); + max_dist = std::max(max_dist, dist); + } + return (max_dist - min_dist) < eps; + } +} // namespace + +VofFlow::DomainInfo::DomainInfo(vtkRectilinearGrid* grid, vtkMPIController* mpiController) + : localExtent_{}, + localZeroExtent_{}, + globalExtent_{}, + localBounds_{}, + localZeroBounds_{}, + globalBounds_{}, + dims_{}, + isUniform_(false), + isParallel_(false), + ghostArray_(nullptr) { + if (grid == nullptr) { + return; + } + + // Save a permanent copy of the grid structure without any data to have access to coords. + grid_ = vtkSmartPointer::New(); + grid_->CopyStructure(grid); + + auto numberOfPieces = grid->GetInformation()->Get(vtkDataObject::DATA_NUMBER_OF_PIECES()); + isParallel_ = mpiController != nullptr && mpiController->GetCommunicator() != nullptr && numberOfPieces > 1; + + // Basic input domain ranges + grid->GetExtent(localExtent_.data()); + grid->GetBounds(localBounds_.data()); + + // cell dimensions (grid->GetDimensions() will return node dimensions) + dims_ = Grid::extentDimensions(localExtent_); + + // cell dimensions + coordsX_ = getArrayValues(grid_->GetXCoordinates()); + coordsY_ = getArrayValues(grid_->GetYCoordinates()); + coordsZ_ = getArrayValues(grid_->GetZCoordinates()); + cellSizesX_ = calcCellSizes(coordsX_); + cellSizesY_ = calcCellSizes(coordsY_); + cellSizesZ_ = calcCellSizes(coordsZ_); + cellCentersX_ = calcCellCenters(coordsX_); + cellCentersY_ = calcCellCenters(coordsY_); + cellCentersZ_ = calcCellCenters(coordsZ_); + + // Validate assumptions for fast ComputeStructuredCoordinates + if (coordsX_.size() < 2 || coordsY_.size() < 2 || coordsZ_.size() < 2) { + throw std::runtime_error("All coordinates arrays must have at least size 2!"); + } + if (!std::is_sorted(coordsX_.begin(), coordsX_.end()) || !std::is_sorted(coordsY_.begin(), coordsY_.end()) || + !std::is_sorted(coordsZ_.begin(), coordsZ_.end())) { + throw std::runtime_error("All coordinates arrays must be sorted!"); + } + + // Uniform Grid + isUniform_ = isUniformCoords(coordsX_) && isUniformCoords(coordsY_) && isUniformCoords(coordsZ_); + + // No MPI / Single process case + // + // Global and zero domain are the same as local domain (no ghost cells and no other threads). + if (!isParallel_) { + globalExtent_ = localExtent_; + localZeroExtent_ = localExtent_; + globalBounds_ = localBounds_; + localZeroBounds_ = localBounds_; + return; + } + + // MPI case extent + // + // GetExtent will include ghost cells, but we also need the local extent without ghost cells ("zero extent"). + // The only official way to find ghost cells seems to lookup within the 'vtkGhostType' array. But here, we want + // to avoid the array traversal and parsing, and calculate the zero extent directly. Unfortunately, only using the + // local + global extent and the number of ghost cells, we cannot reconstruct the zero extent. There are many + // edge cases, where calculations with only these values will fail, i.e. domain boundary reduces number of ghost + // levels. + // Therefore, we copy the VTK internal extent calculation which is used to set up the 'vtkGhostType' array. + // Source: + // https://gitlab.kitware.com/vtk/vtk/-/blob/fc5d7cb1748b7b123c7d97bce88cf9385b39eeb6/Common/ExecutionModel/vtkStreamingDemandDrivenPipeline.cxx#L976-986 + // (VTK commit hash referenced by ParaView v5.9.0 tag) + + grid->GetInformation()->Get(vtkDataObject::ALL_PIECES_EXTENT(), globalExtent_.data()); + auto pieceNumber = grid->GetInformation()->Get(vtkDataObject::DATA_PIECE_NUMBER()); + + vtkExtentTranslator* et = vtkExtentTranslator::New(); + et->PieceToExtentThreadSafe(pieceNumber, numberOfPieces, 0, globalExtent_.data(), localZeroExtent_.data(), + vtkExtentTranslator::BLOCK_MODE, 0); + et->Delete(); + + // MPI case bounds + + vtkDataArray* coordsX = grid->GetXCoordinates(); + vtkDataArray* coordsY = grid->GetYCoordinates(); + vtkDataArray* coordsZ = grid->GetZCoordinates(); + localZeroBounds_[0] = coordsX->GetComponent(localZeroExtent_[0] - localExtent_[0], 0); + localZeroBounds_[1] = coordsX->GetComponent(localZeroExtent_[1] - localExtent_[0], 0); + localZeroBounds_[2] = coordsY->GetComponent(localZeroExtent_[2] - localExtent_[2], 0); + localZeroBounds_[3] = coordsY->GetComponent(localZeroExtent_[3] - localExtent_[2], 0); + localZeroBounds_[4] = coordsZ->GetComponent(localZeroExtent_[4] - localExtent_[4], 0); + localZeroBounds_[5] = coordsZ->GetComponent(localZeroExtent_[5] - localExtent_[4], 0); + + // Save inverse of upper bounds to use combined minimum operation + std::array sendBounds{ + localBounds_[0], + -localBounds_[1], + localBounds_[2], + -localBounds_[3], + localBounds_[4], + -localBounds_[5], + }; + std::array recvBounds{}; + + mpiController->AllReduce(sendBounds.data(), recvBounds.data(), 6, vtkCommunicator::MIN_OP); + + globalBounds_[0] = recvBounds[0]; + globalBounds_[1] = -recvBounds[1]; + globalBounds_[2] = recvBounds[2]; + globalBounds_[3] = -recvBounds[3]; + globalBounds_[4] = recvBounds[4]; + globalBounds_[5] = -recvBounds[5]; + + // Check if isUniform is true over all segments. + // Assuming ghost cells, we can also be sure that there is no slicing change at a process border. + unsigned char sendIsUniform = isUniform_ ? 1 : 0; + unsigned char recvIsUniform = 0; + mpiController->AllReduce(&sendIsUniform, &recvIsUniform, 1, vtkCommunicator::MIN_OP); + isUniform_ = recvIsUniform == 1; + + // While we above avoided the array traversal, the ghost array is still useful to check individual cells without + // calculating if they are within the zero extent. + ghostArray_ = vtkSmartPointer::New(); + ghostArray_->DeepCopy(grid->GetCellGhostArray()); + + // Find neighbors + const int processId = mpiController->GetLocalProcessId(); + const int numProcesses = mpiController->GetNumberOfProcesses(); + + // Send zero extents + std::vector allZeroExtents(6 * numProcesses); + std::vector allZeroBounds(6 * numProcesses); + VofFlow::AllGatherVSameLength(mpiController, localZeroExtent_.data(), allZeroExtents.data(), 6); + VofFlow::AllGatherVSameLength(mpiController, localZeroBounds_.data(), allZeroBounds.data(), 6); + + // Find neighbours + neighborsByDirection_ = std::array, 6>(); + for (int i = 0; i < numProcesses; ++i) { + if (i == processId) { + continue; + } + + extent_t const& zeroExtent_i = *reinterpret_cast(&allZeroExtents[6 * i]); + bounds_t const& zeroBounds_i = *reinterpret_cast(&allZeroBounds[6 * i]); + + // Neighbors are all processes which zero extent overlaps with the ghost cells of the current process. + // We can assume, that zero extents do not overlap. Therefore, we can just intersect the full extent of the + // current process with the zero extent of the other processes. + if (Grid::intersect(localExtent_, zeroExtent_i)) { + // Classify location of each neighbor. + // Store a flag for each of the 6 sides, if the neighbor can be found in this direction. + std::bitset<6> flags; + flags.set(0, zeroExtent_i[0] < localZeroExtent_[0]); // low X + flags.set(1, zeroExtent_i[1] > localZeroExtent_[1]); // up X + flags.set(2, zeroExtent_i[2] < localZeroExtent_[2]); // low Y + flags.set(3, zeroExtent_i[3] > localZeroExtent_[3]); // up Y + flags.set(4, zeroExtent_i[4] < localZeroExtent_[4]); // low Z + flags.set(5, zeroExtent_i[5] > localZeroExtent_[5]); // up Z + + neighbors_.push_back({i, flags, zeroExtent_i, zeroBounds_i}); + + // Convert to by-direction list format + for (int j = 0; j < 6; j++) { + if (flags[j]) { + neighborsByDirection_[j].push_back(i); + } + } + } + } + + // Check if neighbors are symmetric between all processes. This should theoretically always be true, but do some + // paranoid extra error checking here. + // Build a matrix filled with 0 and 1 which processes are neighbors. + std::vector row(numProcesses, 0); + for (const auto& n : neighbors_) { + row[n.processId] = 1; + } + std::vector neighborMatrix(numProcesses * numProcesses, 0); + VofFlow::AllGatherVSameLength(mpiController, row, neighborMatrix); + for (int i = 0; i < numProcesses; i++) { + if (neighborMatrix[i * numProcesses + i] != 0) { + throw std::runtime_error("Bad neighbor finding!"); + } + for (int j = i + 1; j < numProcesses; j++) { + if (neighborMatrix[i * numProcesses + j] != neighborMatrix[j * numProcesses + i]) { + throw std::runtime_error("Bad neighbor finding!"); + } + } + } +} diff --git a/src/VofDataUtils/Grid/DomainInfo.h b/src/VofDataUtils/Grid/DomainInfo.h new file mode 100644 index 0000000..0d2dd42 --- /dev/null +++ b/src/VofDataUtils/Grid/DomainInfo.h @@ -0,0 +1,432 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "../Math/Vector.h" +#include "GridTypes.h" + +class vtkMPIController; + +namespace VofFlow { + class DomainInfo { + public: + struct NeighborInfo { + int processId = 0; + std::bitset<6> flags; + extent_t zeroExtent{}; + bounds_t zeroBounds{}; + }; + + explicit DomainInfo(vtkRectilinearGrid* grid, vtkMPIController* mpiController = nullptr); + + [[nodiscard]] inline vtkRectilinearGrid* grid() const { + return grid_; + } + + [[nodiscard]] inline const extent_t& localExtent() const { + return localExtent_; + } + + [[nodiscard]] inline const extent_t& localZeroExtent() const { + return localZeroExtent_; + } + + [[nodiscard]] inline const extent_t& globalExtent() const { + return globalExtent_; + } + + [[nodiscard]] inline const bounds_t& localBounds() const { + return localBounds_; + } + + [[nodiscard]] inline const bounds_t& localZeroBounds() const { + return localZeroBounds_; + } + + [[nodiscard]] inline const bounds_t& globalBounds() const { + return globalBounds_; + } + + [[nodiscard]] inline const dim_t& cellDims() const { + return dims_; + } + + [[nodiscard]] inline bool isUniform() const { + return isUniform_; + } + + [[nodiscard]] inline bool isParallel() const { + return isParallel_; + } + + [[nodiscard]] inline const vtkSmartPointer& ghostArray() const { + return ghostArray_; + } + + [[nodiscard]] inline const std::vector& coordsX() const { + return coordsX_; + } + + [[nodiscard]] inline const std::vector& coordsY() const { + return coordsY_; + } + + [[nodiscard]] inline const std::vector& coordsZ() const { + return coordsZ_; + } + + [[nodiscard]] inline const std::vector& cellSizesX() const { + return cellSizesX_; + } + + [[nodiscard]] inline const std::vector& cellSizesY() const { + return cellSizesY_; + } + + [[nodiscard]] inline const std::vector& cellSizesZ() const { + return cellSizesZ_; + } + + [[nodiscard]] inline const std::vector& cellCentersX() const { + return cellCentersX_; + } + + [[nodiscard]] inline const std::vector& cellCentersY() const { + return cellCentersY_; + } + + [[nodiscard]] inline const std::vector& cellCentersZ() const { + return cellCentersZ_; + } + + [[nodiscard]] inline const std::vector& neighbors() const { + return neighbors_; + } + + [[nodiscard]] inline const std::array, 6>& neighborsByDirection() const { + return neighborsByDirection_; + } + + [[nodiscard]] inline vtkIdType numCells() const { + return static_cast(dims_[0]) * static_cast(dims_[1]) * + static_cast(dims_[2]); + } + + [[nodiscard]] inline posCoords_t coords(int x, int y, int z) const { + return {coordsX_[x], coordsY_[y], coordsZ_[z]}; + } + + [[nodiscard]] inline posCoords_t coords(gridCoords_t g_coords) const { + return coords(g_coords[0], g_coords[1], g_coords[2]); + } + + [[nodiscard]] inline vec3 coordsVec3(int x, int y, int z) const { + return { + static_cast(coordsX_[x]), + static_cast(coordsY_[y]), + static_cast(coordsZ_[z]), + }; + } + + [[nodiscard]] inline vec3 coordsVec3(gridCoords_t g_coords) const { + return coordsVec3(g_coords[0], g_coords[1], g_coords[2]); + } + + [[nodiscard]] inline posCoords_t cellSize(int g_x, int g_y, int g_z) const { + return {cellSizesX_[g_x], cellSizesY_[g_y], cellSizesZ_[g_z]}; + } + + [[nodiscard]] inline posCoords_t cellSize(gridCoords_t g_coords) const { + return cellSize(g_coords[0], g_coords[1], g_coords[2]); + } + + [[nodiscard]] inline vec3 cellSizeVec3(int g_x, int g_y, int g_z) const { + return { + static_cast(cellSizesX_[g_x]), + static_cast(cellSizesY_[g_y]), + static_cast(cellSizesZ_[g_z]), + }; + } + + [[nodiscard]] inline vec3 cellSizeVec3(gridCoords_t g_coords) const { + return cellSizeVec3(g_coords[0], g_coords[1], g_coords[2]); + } + + [[nodiscard]] inline posCoords_t cellCenter(int g_x, int g_y, int g_z) const { + return {cellCentersX_[g_x], cellCentersY_[g_y], cellCentersZ_[g_z]}; + } + + [[nodiscard]] inline posCoords_t cellCenter(gridCoords_t g_coords) const { + return cellCenter(g_coords[0], g_coords[1], g_coords[2]); + } + + [[nodiscard]] inline vec3 cellCenterVec3(int g_x, int g_y, int g_z) const { + return { + static_cast(cellCentersX_[g_x]), + static_cast(cellCentersY_[g_y]), + static_cast(cellCentersZ_[g_z]), + }; + } + + [[nodiscard]] inline vec3 cellCenterVec3(gridCoords_t g_coords) const { + return cellCenterVec3(g_coords[0], g_coords[1], g_coords[2]); + } + + [[nodiscard]] inline bool isInLocalExtent(const gridCoords_t& e_coords) const { + return Grid::isInExtent(localExtent_, e_coords); + } + + [[nodiscard]] inline bool isInLocalExtent(int e_x, int e_y, int e_z) const { + return isInLocalExtent(gridCoords_t{e_x, e_y, e_z}); + } + + [[nodiscard]] inline bool isInZeroExtent(const gridCoords_t& e_coords) const { + return Grid::isInExtent(localZeroExtent_, e_coords); + } + + [[nodiscard]] inline bool isInZeroExtent(int e_x, int e_y, int e_z) const { + return isInZeroExtent(gridCoords_t{e_x, e_y, e_z}); + } + + [[nodiscard]] inline bool isInZeroExtent(vtkIdType idx) const { + return isInZeroExtent(idxToExtentCoord(idx)); + } + + [[nodiscard]] inline vtkIdType gridCoordToIdx(int g_x, int g_y, int g_z) const { + return static_cast(g_x) + static_cast(g_y) * static_cast(dims_[0]) + + static_cast(g_z) * static_cast(dims_[0]) * static_cast(dims_[1]); + } + + [[nodiscard]] inline vtkIdType gridCoordToIdx(const gridCoords_t& g_coords) const { + return gridCoordToIdx(g_coords[0], g_coords[1], g_coords[2]); + } + + [[nodiscard]] inline gridCoords_t idxToGridCoord(vtkIdType idx) const { + const int g_x = static_cast(idx % dims_[0]); + const int g_y = static_cast((idx / dims_[0]) % dims_[1]); + const int g_z = static_cast(idx / (dims_[0] * dims_[1])); + return {g_x, g_y, g_z}; + } + + [[nodiscard]] inline vtkIdType localExtentCoordToIdx(int e_x, int e_y, int e_z) const { + if (!isInLocalExtent(e_x, e_y, e_z)) { + throw std::runtime_error("Extent coord out of bounds!"); + } + return (e_x - localExtent_[0]) + (e_y - localExtent_[2]) * dims_[0] + + (e_z - localExtent_[4]) * dims_[0] * dims_[1]; + } + + [[nodiscard]] inline vtkIdType localExtentCoordToIdx(const gridCoords_t& e_coords) const { + return localExtentCoordToIdx(e_coords[0], e_coords[1], e_coords[2]); + } + + [[nodiscard]] inline gridCoords_t idxToExtentCoord(vtkIdType idx) const { + const auto [g_x, g_y, g_z] = idxToGridCoord(idx); + return {localExtent_[0] + g_x, localExtent_[2] + g_y, localExtent_[4] + g_z}; + } + + [[nodiscard]] inline gridCoords_t localExtentCoordToGridCoord(int e_x, int e_y, int e_z) const { + if (!isInLocalExtent(e_x, e_y, e_z)) { + throw std::runtime_error("Extent coord out of bounds!"); + } + return {e_x - localExtent_[0], e_y - localExtent_[2], e_z - localExtent_[4]}; + } + + [[nodiscard]] inline gridCoords_t localExtentCoordToGridCoord(const gridCoords_t& e_coords) const { + return localExtentCoordToGridCoord(e_coords[0], e_coords[1], e_coords[2]); + } + + [[nodiscard]] inline bool isInLocalBounds(const posCoords_t& pos) const { + return Grid::isInBounds(localBounds_, pos); + } + + [[nodiscard]] inline bool isInLocalBounds(double x, double y, double z) const { + return isInLocalBounds(posCoords_t{x, y, z}); + } + + [[nodiscard]] inline bool isInLocalBounds(const vec3& p) const { + return isInLocalBounds(p.x, p.y, p.z); + } + + [[nodiscard]] inline bool isInZeroBounds(const posCoords_t& pos) const { + return Grid::isInBounds(localZeroBounds_, pos); + } + + [[nodiscard]] inline bool isInZeroBounds(double x, double y, double z) const { + return isInZeroBounds(posCoords_t{x, y, z}); + } + + [[nodiscard]] inline bool isInZeroBounds(const vec3& p) const { + return isInZeroBounds(p.x, p.y, p.z); + } + + [[nodiscard]] inline bool isInGlobalBounds(const posCoords_t& pos) const { + return Grid::isInBounds(globalBounds_, pos); + } + + [[nodiscard]] inline bool isInGlobalBounds(double x, double y, double z) const { + return isInGlobalBounds(posCoords_t{x, y, z}); + } + + [[nodiscard]] inline bool isInGlobalBounds(const vec3& p) const { + return isInGlobalBounds(p.x, p.y, p.z); + } + + [[nodiscard]] inline int isOutOfZeroBoundsDirection(const vec3& p) const { + return Grid::isOutOfBoundsDirection(localZeroBounds_, posCoords_t{p.x, p.y, p.z}); + } + + [[nodiscard]] inline std::optional computeStructuredCoordinates( + const posCoords_t& p) const { + StructuredCoordinates sc; + if (computeStructuredCoordinates(coordsX_, p[0], sc.cellCoords[0], sc.relativeCoords[0]) && + computeStructuredCoordinates(coordsY_, p[1], sc.cellCoords[1], sc.relativeCoords[1]) && + computeStructuredCoordinates(coordsZ_, p[2], sc.cellCoords[2], sc.relativeCoords[2])) { + return sc; + } + return std::nullopt; + } + + [[nodiscard]] inline std::optional computeStructuredCoordinates(const vec3& p) const { + return computeStructuredCoordinates(posCoords_t{p.x, p.y, p.z}); + } + + [[nodiscard]] inline StructuredCoordinates computeNearestStructuredCoordinates(const posCoords_t& p) const { + StructuredCoordinates sc; + computeStructuredCoordinates(coordsX_, p[0], sc.cellCoords[0], sc.relativeCoords[0]); + computeStructuredCoordinates(coordsY_, p[1], sc.cellCoords[1], sc.relativeCoords[1]); + computeStructuredCoordinates(coordsZ_, p[2], sc.cellCoords[2], sc.relativeCoords[2]); + return sc; + } + + [[nodiscard]] inline StructuredCoordinates computeNearestStructuredCoordinates(const vec3& p) const { + return computeNearestStructuredCoordinates(posCoords_t{p.x, p.y, p.z}); + } + + [[nodiscard]] inline std::optional posToGridCoord(const posCoords_t& p) const { + auto coords = computeStructuredCoordinates(p); + if (coords.has_value()) { + return coords.value().cellCoords; + } + return std::nullopt; + } + + [[nodiscard]] inline std::optional posToGridCoord(const vec3& p) const { + return posToGridCoord(posCoords_t{p.x, p.y, p.z}); + } + + [[nodiscard]] inline gridCoords_t posToGridCoordEx(const posCoords_t& p) const { + const auto cell = posToGridCoord(p); + if (cell.has_value()) { + return cell.value(); + } + throw std::runtime_error("Position is outside of grid bounds!"); + } + + [[nodiscard]] inline gridCoords_t posToGridCoordEx(const vec3& p) const { + return posToGridCoordEx(posCoords_t{p.x, p.y, p.z}); + } + + [[nodiscard]] inline std::optional posToIdx(const posCoords_t& p) const { + const auto g_coords = posToGridCoord(p); + if (g_coords.has_value()) { + return gridCoordToIdx(g_coords.value()); + } + return std::nullopt; + } + + [[nodiscard]] inline std::optional posToIdx(const vec3& p) const { + return posToIdx(posCoords_t{p.x, p.y, p.z}); + } + + [[nodiscard]] inline vtkIdType posToIdxEx(const posCoords_t& p) const { + return gridCoordToIdx(posToGridCoordEx(p)); + } + + [[nodiscard]] inline vtkIdType posToIdxEx(const vec3& p) const { + return gridCoordToIdx(posToGridCoordEx(p)); + } + + private: + /** + * VTK's `ComputeStructuredCoordinates()` is very slow. The two major points are the linear search and the + * overhead from the GetComponent interface of the arrays. Here, we implement binary search and directly use + * the std::vector coords. + * Assumptions: coords is sorted and size >= 2. (Both checked in constructor.) + */ + inline bool computeStructuredCoordinates(const std::vector& coords, double pos, int& cell, + double& pcoord) const { + // If out of bounds return closest point. + if (pos < coords.front()) { + cell = 0; + pcoord = 0.0; + return false; + } else if (pos >= coords.back()) { + cell = static_cast(coords.size()) - 2; + pcoord = 1.0; + return false; + } + + // Shortcut for uniform grids + if (isUniform_) { + double relative_pos = (pos - coords.front()) / (coords.back() - coords.front()); // [0, 1) + relative_pos *= static_cast(coords.size() - 1); // [0, num_cells) + double iptr = 0.0; + pcoord = std::clamp(std::modf(relative_pos, &iptr), 0.0, 1.0); + cell = std::clamp(static_cast(iptr), 0, static_cast(coords.size() - 1)); + return true; + } + + std::size_t left = 0; + std::size_t right = coords.size() - 1; + while (left + 1 < right) { + const std::size_t mid = left + (right - left) / 2; + if (pos < coords[mid]) { + right = mid; + } else { + left = mid; + } + } + cell = static_cast(left); + pcoord = (pos - coords[left]) / (coords[right] - coords[left]); + return true; + } + + vtkSmartPointer grid_; + extent_t localExtent_; // Extent of thread local grid. + extent_t localZeroExtent_; // Extent of the thread local grid without ghost cells. + extent_t globalExtent_; // Extent of the global grid across all threads. + bounds_t localBounds_; // Bounds of the thread local grid domain. + bounds_t localZeroBounds_; // Bounds of the thread local grid domain without ghost cells. + bounds_t globalBounds_; // Bounds of the global grid domain across all threads. + dim_t dims_; + bool isUniform_; + bool isParallel_; + vtkSmartPointer ghostArray_; + + std::vector coordsX_; + std::vector coordsY_; + std::vector coordsZ_; + std::vector cellSizesX_; + std::vector cellSizesY_; + std::vector cellSizesZ_; + std::vector cellCentersX_; + std::vector cellCentersY_; + std::vector cellCentersZ_; + + std::vector neighbors_; + std::array, 6> neighborsByDirection_; // list neighbors by direction + }; +} // namespace VofFlow diff --git a/src/VofDataUtils/Grid/Gradient.cpp b/src/VofDataUtils/Grid/Gradient.cpp new file mode 100644 index 0000000..2634af0 --- /dev/null +++ b/src/VofDataUtils/Grid/Gradient.cpp @@ -0,0 +1,199 @@ +#include "Gradient.h" + +#include + +#include +#include +#include +#include + +namespace { + struct GradientWorker { + const VofFlow::DomainInfo& domainInfo_; + const int x_; + const int y_; + const int z_; + VofFlow::vec3 gradient_; + const bool borderXmin_; + const bool borderXmax_; + const bool borderYmin_; + const bool borderYmax_; + const bool borderZmin_; + const bool borderZmax_; + + GradientWorker(const VofFlow::DomainInfo& domainInfo, const VofFlow::gridCoords_t& g_coords) + : domainInfo_(domainInfo), + x_(g_coords[0]), + y_(g_coords[1]), + z_(g_coords[2]), + gradient_(0.0f, 0.0f, 0.0f), + borderXmin_(x_ <= 0), + borderXmax_(x_ >= domainInfo_.cellDims()[0] - 1), + borderYmin_(y_ <= 0), + borderYmax_(y_ >= domainInfo_.cellDims()[1] - 1), + borderZmin_(z_ <= 0), + borderZmax_(z_ >= domainInfo_.cellDims()[2] - 1) {} + + template + inline ValueType gradientXDim(const vtkDataArrayAccessor>& data, int y, + int z) { + const auto& cellSizesX = domainInfo_.cellSizesX(); + if (borderXmin_) { + // forward difference + const ValueType val0 = data.Get(domainInfo_.gridCoordToIdx(x_, y, z), 0); + const ValueType val1 = data.Get(domainInfo_.gridCoordToIdx(x_ + 1, y, z), 0); + return (val1 - val0) / (0.5 * (cellSizesX[x_] + cellSizesX[x_ + 1])); + } else if (borderXmax_) { + // backward difference + const ValueType val0 = data.Get(domainInfo_.gridCoordToIdx(x_ - 1, y, z), 0); + const ValueType val1 = data.Get(domainInfo_.gridCoordToIdx(x_, y, z), 0); + return (val1 - val0) / (0.5 * (cellSizesX[x_ - 1] + cellSizesX[x_])); + } else { + // central difference + const ValueType val0 = data.Get(domainInfo_.gridCoordToIdx(x_ - 1, y, z), 0); + const ValueType val1 = data.Get(domainInfo_.gridCoordToIdx(x_ + 1, y, z), 0); + return (val1 - val0) / (0.5 * cellSizesX[x_ - 1] + cellSizesX[x_] + 0.5 * cellSizesX[x_ + 1]); + } + } + + template + inline ValueType gradientYDim(const vtkDataArrayAccessor>& data, int x, + int z) { + const auto& cellSizesY = domainInfo_.cellSizesY(); + if (borderYmin_) { + // forward difference + const ValueType val0 = data.Get(domainInfo_.gridCoordToIdx(x, y_, z), 0); + const ValueType val1 = data.Get(domainInfo_.gridCoordToIdx(x, y_ + 1, z), 0); + return (val1 - val0) / (0.5 * (cellSizesY[y_] + cellSizesY[y_ + 1])); + } else if (borderYmax_) { + // backward difference + const ValueType val0 = data.Get(domainInfo_.gridCoordToIdx(x, y_ - 1, z), 0); + const ValueType val1 = data.Get(domainInfo_.gridCoordToIdx(x, y_, z), 0); + return (val1 - val0) / (0.5 * (cellSizesY[y_ - 1] + cellSizesY[y_])); + } else { + // central difference + const ValueType val0 = data.Get(domainInfo_.gridCoordToIdx(x, y_ - 1, z), 0); + const ValueType val1 = data.Get(domainInfo_.gridCoordToIdx(x, y_ + 1, z), 0); + return (val1 - val0) / (0.5 * cellSizesY[y_ - 1] + cellSizesY[y_] + 0.5 * cellSizesY[y_ + 1]); + } + } + + template + inline ValueType gradientZDim(const vtkDataArrayAccessor>& data, int x, + int y) { + const auto& cellSizesZ = domainInfo_.cellSizesZ(); + if (borderZmin_) { + // forward difference + const ValueType val0 = data.Get(domainInfo_.gridCoordToIdx(x, y, z_), 0); + const ValueType val1 = data.Get(domainInfo_.gridCoordToIdx(x, y, z_ + 1), 0); + return (val1 - val0) / (0.5 * (cellSizesZ[z_] + cellSizesZ[z_ + 1])); + } else if (borderZmax_) { + // backward difference + const ValueType val0 = data.Get(domainInfo_.gridCoordToIdx(x, y, z_ - 1), 0); + const ValueType val1 = data.Get(domainInfo_.gridCoordToIdx(x, y, z_), 0); + return (val1 - val0) / (0.5 * (cellSizesZ[z_ - 1] + cellSizesZ[z_])); + } else { + // central difference + const ValueType val0 = data.Get(domainInfo_.gridCoordToIdx(x, y, z_ - 1), 0); + const ValueType val1 = data.Get(domainInfo_.gridCoordToIdx(x, y, z_ + 1), 0); + return (val1 - val0) / (0.5 * cellSizesZ[z_ - 1] + cellSizesZ[z_] + 0.5 * cellSizesZ[z_ + 1]); + } + } + + template + void operator()(vtkAOSDataArrayTemplate* dataArray) { + VTK_ASSUME(dataArray->GetNumberOfComponents() == 1); + + vtkDataArrayAccessor> data(dataArray); + + // TODO weighting by cell size is not implemented yet! + + // X + ValueType gradientX = 0.0; + { + ValueType weightsX = 0.0; + for (int y = 0; y < 3; y++) { + if ((y == 0 && borderYmin_) || (y == 2 && borderYmax_)) { + continue; + } + const int weightY = (y == 1) ? 2 : 1; + for (int z = 0; z < 3; z++) { + if ((z == 0 && borderZmin_) || (z == 2 && borderZmax_)) { + continue; + } + const int weightZ = (z == 1) ? 2 : 1; + const double weight = static_cast(weightY * weightZ); + + gradientX += weight * gradientXDim(data, y_ - 1 + y, z_ - 1 + z); + weightsX += weight; + } + } + gradientX /= weightsX; + } + + // Y + ValueType gradientY = 0.0; + { + ValueType weightsY = 0.0; + for (int x = 0; x < 3; x++) { + if ((x == 0 && borderXmin_) || (x == 2 && borderXmax_)) { + continue; + } + const int weightX = (x == 1) ? 2 : 1; + for (int z = 0; z < 3; z++) { + if ((z == 0 && borderZmin_) || (z == 2 && borderZmax_)) { + continue; + } + const int weightZ = (z == 1) ? 2 : 1; + const double weight = static_cast(weightX * weightZ); + + gradientY += weight * gradientYDim(data, x_ - 1 + x, z_ - 1 + z); + weightsY += weight; + } + } + gradientY /= weightsY; + } + + // Z + ValueType gradientZ = 0.0; + { + ValueType weightsZ = 0.0; + for (int x = 0; x < 3; x++) { + if ((x == 0 && borderXmin_) || (x == 2 && borderXmax_)) { + continue; + } + const int weightX = (x == 1) ? 2 : 1; + for (int y = 0; y < 3; y++) { + if ((y == 0 && borderYmin_) || (y == 2 && borderYmax_)) { + continue; + } + const int weightY = (y == 1) ? 2 : 1; + const double weight = static_cast(weightX * weightY); + + gradientZ += weight * gradientZDim(data, x_ - 1 + x, y_ - 1 + y); + weightsZ += weight; + } + } + gradientZ /= weightsZ; + } + + gradient_ = VofFlow::vec3{ + static_cast(gradientX), + static_cast(gradientY), + static_cast(gradientZ), + }; + } + }; +} // namespace + +VofFlow::vec3 VofFlow::gradient(const DomainInfo& domainInfo, const gridCoords_t& g_coords, + const vtkSmartPointer& data) { + GradientWorker worker(domainInfo, g_coords); + + typedef vtkArrayDispatch::DispatchByValueType Dispatcher; + if (!Dispatcher::Execute(data, worker)) { + throw std::runtime_error("Cannot dispatch array worker!"); + } + + return worker.gradient_; +} diff --git a/src/VofDataUtils/Grid/Gradient.h b/src/VofDataUtils/Grid/Gradient.h new file mode 100644 index 0000000..2a50f6c --- /dev/null +++ b/src/VofDataUtils/Grid/Gradient.h @@ -0,0 +1,12 @@ +#pragma once + +#include +#include + +#include "../Math/Vector.h" +#include "DomainInfo.h" +#include "GridTypes.h" + +namespace VofFlow { + vec3 gradient(const DomainInfo& domain, const gridCoords_t& g_coords, const vtkSmartPointer& data); +} diff --git a/src/VofDataUtils/Grid/GridData.h b/src/VofDataUtils/Grid/GridData.h new file mode 100644 index 0000000..6bb2a77 --- /dev/null +++ b/src/VofDataUtils/Grid/GridData.h @@ -0,0 +1,44 @@ +#pragma once + +#include +#include + +#include +#include + +#include "../VtkUtil/VtkUtilArray.h" +#include "DomainInfo.h" +#include "GridTypes.h" + +namespace VofFlow { + template + vtkSmartPointer extractExtent(const DomainInfo& domainInfo, extent_t extent, + vtkSmartPointer array) { + // Validate requested extent is fully within the array extent. + if (extent[0] < domainInfo.localExtent()[0] || extent[1] > domainInfo.localExtent()[1] || + extent[2] < domainInfo.localExtent()[2] || extent[3] > domainInfo.localExtent()[3] || + extent[4] < domainInfo.localExtent()[4] || extent[5] > domainInfo.localExtent()[5]) { + throw std::runtime_error("Requested extent is out of bounds!"); + } + + const int numComponents = array->GetNumberOfComponents(); + const auto extDim = Grid::extentDimensions(extent); + + vtkSmartPointer result = vtkSmartPointer::New(); + result->SetName(array->GetName()); + result->SetNumberOfComponents(numComponents); + result->SetNumberOfTuples(Grid::extentNumCells(extent)); + + for (int z = 0; z < extDim[2]; z++) { + for (int y = 0; y < extDim[1]; y++) { + vtkIdType fromIdx = domainInfo.localExtentCoordToIdx(extent[0], extent[2] + y, extent[4] + z); + vtkIdType toIdx = z * extDim[1] * extDim[0] + y * extDim[0]; + auto* fromPtr = getTuplePointer(array, fromIdx); + auto* toPtr = getTuplePointer(result, toIdx); + std::memcpy(toPtr, fromPtr, extDim[0] * numComponents * sizeof(typename ArrayT::ValueType)); + } + } + + return result; + } +} // namespace VofFlow diff --git a/src/VofDataUtils/Grid/GridIterator.h b/src/VofDataUtils/Grid/GridIterator.h new file mode 100644 index 0000000..78101b7 --- /dev/null +++ b/src/VofDataUtils/Grid/GridIterator.h @@ -0,0 +1,88 @@ +#pragma once + +#include + +#include "GridTypes.h" + +namespace VofFlow { + class GridRange { + public: + GridRange(const gridCoords_t& startPos, const gridCoords_t& endPos) : startPos_(startPos), endPos_(endPos) {} + + explicit GridRange(const extent_t& extent) { + startPos_ = {extent[0], extent[2], extent[4]}; + endPos_ = {extent[1], extent[3], extent[5]}; + } + + GridRange(const dim_t& dimensions, const gridCoords_t& start, int searchRange) { + startPos_ = { + std::max(start[0] - searchRange, 0), + std::max(start[1] - searchRange, 0), + std::max(start[2] - searchRange, 0), + }; + endPos_ = { + std::min(start[0] + searchRange + 1, dimensions[0]), + std::min(start[1] + searchRange + 1, dimensions[1]), + std::min(start[2] + searchRange + 1, dimensions[2]), + }; + } + + class Iterator { + public: + Iterator(const gridCoords_t& startPos, const gridCoords_t& endPos) + : position_(startPos), + startPos_(startPos), + endPos_(endPos) { + if (position_[0] >= endPos_[0] || position_[1] >= endPos_[1] || position_[2] >= endPos_[2]) { + position_ = endPos_; + } + } + + inline Iterator& operator++() { + position_[0]++; + if (position_[0] >= endPos_[0]) { + position_[0] = startPos_[0]; + position_[1]++; + if (position_[1] >= endPos_[1]) { + position_[1] = startPos_[1]; + position_[2]++; + if (position_[2] >= endPos_[2]) { + position_ = endPos_; + } + } + } + return *this; + } + + inline const gridCoords_t& operator*() const { + return position_; + } + + inline bool operator==(const Iterator& other) { + return position_ == other.position_; + } + + inline bool operator!=(const Iterator& other) { + return position_ != other.position_; + } + + private: + gridCoords_t position_; + const gridCoords_t& startPos_; + const gridCoords_t& endPos_; + }; + + [[nodiscard]] inline auto begin() const { + return Iterator(startPos_, endPos_); + } + + [[nodiscard]] inline auto end() const { + return Iterator(endPos_, endPos_); + } + + private: + gridCoords_t startPos_; + gridCoords_t endPos_; + }; + +} // namespace VofFlow diff --git a/src/VofDataUtils/Grid/GridTypes.h b/src/VofDataUtils/Grid/GridTypes.h new file mode 100644 index 0000000..7e16b34 --- /dev/null +++ b/src/VofDataUtils/Grid/GridTypes.h @@ -0,0 +1,94 @@ +#pragma once + +#include +#include +#include + +namespace VofFlow { + using extent_t = std::array; + using bounds_t = std::array; + using dim_t = std::array; + using gridCoords_t = std::array; + using posCoords_t = std::array; + + struct StructuredCoordinates { + gridCoords_t cellCoords{}; + posCoords_t relativeCoords{}; + }; + + namespace Grid { + inline bool isValidExtent(const extent_t& e) { + return e[0] < e[1] && e[2] < e[3] && e[4] < e[5]; + } + + inline bool isInExtent(const extent_t& e, const gridCoords_t& c) { + return c[0] >= e[0] && c[0] < e[1] && c[1] >= e[2] && c[1] < e[3] && c[2] >= e[4] && c[2] < e[5]; + } + + inline dim_t extentDimensions(const extent_t& e) { + if (isValidExtent(e)) { + return {e[1] - e[0], e[3] - e[2], e[5] - e[4]}; + } + return {0, 0, 0}; + } + + inline int extentNumCells(const extent_t& e) { + const auto& dim = extentDimensions(e); + return dim[0] * dim[1] * dim[2]; + } + + inline bool intersect(const extent_t& a, const extent_t& b) { + // VTK extent upper border marks first cell after the grid extent. Therefore, use '<' instead of '<=' to + // only detect actually overlapping grid extents and not just 'touching' extents. + return a[0] < b[1] && a[1] > b[0] && a[2] < b[3] && a[3] > b[2] && a[4] < b[5] && a[5] > b[4]; + } + + inline extent_t intersection(const extent_t& a, const extent_t& b) { + return { + std::max(a[0], b[0]), + std::min(a[1], b[1]), + std::max(a[2], b[2]), + std::min(a[3], b[3]), + std::max(a[4], b[4]), + std::min(a[5], b[5]), + }; + } + + inline bool isInBounds(const bounds_t& b, const posCoords_t& p) { + return p[0] >= b[0] && p[0] < b[1] && p[1] >= b[2] && p[1] < b[3] && p[2] >= b[4] && p[2] < b[5]; + } + + inline int isOutOfBoundsDirection(const bounds_t& b, const posCoords_t& p) { + if (p[0] < b[0]) { + return 0; + } + if (p[0] >= b[1]) { + return 1; + } + if (p[1] < b[2]) { + return 2; + } + if (p[1] >= b[3]) { + return 3; + } + if (p[2] < b[4]) { + return 4; + } + if (p[2] >= b[5]) { + return 5; + } + return -1; + } + + template + inline std::string toString(const std::array& a) { + return std::to_string(a[0]) + " " + std::to_string(a[1]) + " " + std::to_string(a[2]) + " " + + std::to_string(a[3]) + " " + std::to_string(a[4]) + " " + std::to_string(a[5]); + } + + template + inline std::string toString(const std::array& a) { + return std::to_string(a[0]) + " " + std::to_string(a[1]) + " " + std::to_string(a[2]); + } + } // namespace Grid +} // namespace VofFlow diff --git a/src/VofDataUtils/Math/Vector.h b/src/VofDataUtils/Math/Vector.h new file mode 100644 index 0000000..025e348 --- /dev/null +++ b/src/VofDataUtils/Math/Vector.h @@ -0,0 +1,192 @@ +#pragma once + +#include +#include +#include + +namespace VofFlow { + struct vec3 { + float x; + float y; + float z; + + vec3() : x(0.0f), y(0.0f), z(0.0f){}; + explicit vec3(float v) : x(v), y(v), z(v){}; + explicit vec3(float v[3]) : x(v[0]), y(v[1]), z(v[2]){}; + vec3(float x, float y, float z) : x(x), y(y), z(z){}; + + // Allow float[3] style pointer access. + [[nodiscard]] inline const float* data() const { + static_assert(sizeof(vec3) == 3 * sizeof(float)); + return reinterpret_cast(this); + } + + inline static vec3 zero() { + return {0.0f, 0.0f, 0.0f}; + } + + inline static vec3 one() { + return {1.0f, 1.0f, 1.0f}; + } + }; + + inline vec3 operator-(const vec3& a) { + return {-a.x, -a.y, -a.z}; + } + + inline vec3 operator+(const vec3& a, const vec3& b) { + return {a.x + b.x, a.y + b.y, a.z + b.z}; + } + + inline void operator+=(vec3& a, const vec3& b) { + a.x += b.x; + a.y += b.y; + a.z += b.z; + } + + inline vec3 operator+(const vec3& a, float b) { + return {a.x + b, a.y + b, a.z + b}; + } + + inline void operator+=(vec3& a, float b) { + a.x += b; + a.y += b; + a.z += b; + } + + inline vec3 operator-(const vec3& a, const vec3& b) { + return {a.x - b.x, a.y - b.y, a.z - b.z}; + } + + inline void operator-=(vec3& a, const vec3& b) { + a.x -= b.x; + a.y -= b.y; + a.z -= b.z; + } + + inline vec3 operator-(const vec3& a, float b) { + return {a.x - b, a.y - b, a.z - b}; + } + + inline void operator-=(vec3& a, float b) { + a.x -= b; + a.y -= b; + a.z -= b; + } + + inline vec3 operator*(const vec3& a, const vec3& b) { + return {a.x * b.x, a.y * b.y, a.z * b.z}; + } + + inline void operator*=(vec3& a, const vec3& b) { + a.x *= b.x; + a.y *= b.y; + a.z *= b.z; + } + + inline vec3 operator*(const vec3& a, float b) { + return {a.x * b, a.y * b, a.z * b}; + } + + inline void operator*=(vec3& a, float b) { + a.x *= b; + a.y *= b; + a.z *= b; + } + + inline vec3 operator*(float a, const vec3& b) { + return {a * b.x, a * b.y, a * b.z}; + } + + inline vec3 operator/(const vec3& a, const vec3& b) { + return {a.x / b.x, a.y / b.y, a.z / b.z}; + } + + inline void operator/=(vec3& a, const vec3& b) { + a.x /= b.x; + a.y /= b.y; + a.z /= b.z; + } + + inline vec3 operator/(const vec3& a, float b) { + return {a.x / b, a.y / b, a.z / b}; + } + + inline void operator/=(vec3& a, float b) { + a.x /= b; + a.y /= b; + a.z /= b; + } + + inline float dot(const vec3& a, const vec3& b) { + return a.x * b.x + a.y * b.y + a.z * b.z; + } + + inline float lengthSq(const vec3& v) { + return dot(v, v); + } + + inline float length(const vec3& v) { + return std::sqrt(lengthSq(v)); + } + + inline vec3 normalize(const vec3& v) { + return v / length(v); + } + + inline float distanceSq(const vec3& a, const vec3& b) { + return lengthSq(b - a); + } + + inline float distance(const vec3& a, const vec3& b) { + return std::sqrt(distanceSq(a, b)); + } + + inline vec3 cross(const vec3& a, const vec3& b) { + return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x}; + } + + inline vec3 lerp(const vec3& a, const vec3& b, float t) { + return a + t * (b - a); + } + + inline vec3 clamp(const vec3& v, const vec3& lo, const vec3& hi) { + return { + std::clamp(v.x, lo.x, hi.x), + std::clamp(v.y, lo.y, hi.y), + std::clamp(v.z, lo.z, hi.z), + }; + } + + inline vec3 min(const vec3& a, const vec3& b) { + return { + std::min(a.x, b.x), + std::min(a.y, b.y), + std::min(a.z, b.z), + }; + } + + inline vec3 max(const vec3& a, const vec3& b) { + return { + std::max(a.x, b.x), + std::max(a.y, b.y), + std::max(a.z, b.z), + }; + } + + inline vec3 nextafter(const vec3& from, const vec3& to) { + return { + std::nextafter(from.x, to.x), + std::nextafter(from.y, to.y), + std::nextafter(from.z, to.z), + }; + } + + inline bool testInRange(const vec3& v, const vec3& lo, const vec3& hi) { + return v.x >= lo.x && v.x <= hi.x && v.y >= lo.y && v.y <= hi.y && v.z >= lo.z && v.z <= hi.z; + } + + inline std::ostream& operator<<(std::ostream& os, const vec3& v) { + return os << "(" << v.x << " " << v.y << " " << v.z << ")"; + } +} // namespace VofFlow diff --git a/src/VofDataUtils/Misc/CgalUtil.h b/src/VofDataUtils/Misc/CgalUtil.h new file mode 100644 index 0000000..d8e16f5 --- /dev/null +++ b/src/VofDataUtils/Misc/CgalUtil.h @@ -0,0 +1,52 @@ +#pragma once + +#include +#include + +#include +#include +#include +#include + +#include "../Grid/GridTypes.h" +#include "../Math/Vector.h" + +namespace VofFlow { + typedef CGAL::Exact_predicates_inexact_constructions_kernel K; + + inline static K::Point_3 toPoint_3(const posCoords_t& p) { + return K::Point_3(p[0], p[1], p[2]); + } + + inline static K::Point_3 toPoint_3(const vec3& p) { + return K::Point_3(p.x, p.y, p.z); + } + + inline static vec3 toVec3(const K::Point_3& p) { + return { + static_cast(CGAL::to_double(p.x())), + static_cast(CGAL::to_double(p.y())), + static_cast(CGAL::to_double(p.z())), + }; + } + + inline static std::vector toVec3(const std::vector& points) { + std::vector result(points.size()); + for (std::size_t i = 0; i < points.size(); i++) { + result[i] = toVec3(points[i]); + } + return result; + } + + inline static K::Vector_3 toVector_3(const vec3& v) { + return K::Vector_3(v.x, v.y, v.z); + } + + inline static vec3 toVec3(const K::Vector_3& v) { + return { + static_cast(CGAL::to_double(v.x())), + static_cast(CGAL::to_double(v.y())), + static_cast(CGAL::to_double(v.z())), + }; + } +} // namespace VofFlow diff --git a/src/VofDataUtils/Misc/ListSearch.h b/src/VofDataUtils/Misc/ListSearch.h new file mode 100644 index 0000000..49c29a3 --- /dev/null +++ b/src/VofDataUtils/Misc/ListSearch.h @@ -0,0 +1,25 @@ +#pragma once + +#include +#include +#include +#include + +namespace VofFlow { + template + std::size_t findClosestIndex(const T value, const std::vector& list) { + if (list.empty()) { + throw std::runtime_error("List is empty!"); + } + std::size_t idx = 0; + T minDist = std::abs(list[0] - value); + for (std::size_t i = 1; i < list.size(); i++) { + const T dist = std::abs(list[i] - value); + if (dist < minDist) { + minDist = dist; + idx = i; + } + } + return idx; + } +} // namespace VofFlow diff --git a/src/VofDataUtils/Misc/Profiling.h b/src/VofDataUtils/Misc/Profiling.h new file mode 100644 index 0000000..ba18ed7 --- /dev/null +++ b/src/VofDataUtils/Misc/Profiling.h @@ -0,0 +1,11 @@ +#pragma once + +#ifdef VOFFLOW_USE_TRACY +#include +#else + +// Use empty macro definitions from Tracy.hpp +#define ZoneScoped +#define ZoneName(x, y) + +#endif diff --git a/src/VofDataUtils/Misc/VofData.h b/src/VofDataUtils/Misc/VofData.h new file mode 100644 index 0000000..68f734d --- /dev/null +++ b/src/VofDataUtils/Misc/VofData.h @@ -0,0 +1,21 @@ +#pragma once + +#include +#include + +namespace VofFlow { + /* + * In two phases/species data there is only a single field usually called `vof` (`funs*.hdf`). In three + * phases/species data there are two fields called `vof` (`funs*.hdf`) and `f3` (`fun3*.hdf`). Here, the PLIC + * surface of the `f3` field is based on the gradient, while the PLIC surface of the `vof` field uses additional + * normals in the dataset `n_c_3ph` (`n3ph*.hdf`). The `f3` field in the three phases/species case and the `vof` + * field in the two phases/species case are handled the same way, despite the name being different. + * Therefore, we use the terms `vof1st` and `vof2nd` to simplify the different cases. So `vof1st` is always the + * field using the gradient, while `vof2nd` (if present) is the grid using the additional normals. + */ + struct VofData { + vtkSmartPointer vof1st; + vtkSmartPointer vof2nd; + vtkSmartPointer vof2ndNorm; + }; +} // namespace VofFlow diff --git a/src/VofDataUtils/Plic/CachedPlic.h b/src/VofDataUtils/Plic/CachedPlic.h new file mode 100644 index 0000000..6e276a2 --- /dev/null +++ b/src/VofDataUtils/Plic/CachedPlic.h @@ -0,0 +1,48 @@ +#pragma once + +#include +#include +#include + +#include + +#include "../Grid/DomainInfo.h" +#include "../Misc/Profiling.h" +#include "../Misc/VofData.h" +#include "PlicUtil.h" + +namespace VofFlow { + + class CachedPlic { + public: + explicit CachedPlic(const DomainInfo& domainInfo, const VofData& vofData, double eps, std::size_t numIterations) + : domainInfo_(domainInfo), + vofData_(vofData), + eps_(eps), + numIterations_(numIterations){}; + + inline const PlicCellClassResult& get(vtkIdType gridIdx) { + ZoneScoped; + + auto it = cache_.find(gridIdx); + if (it == cache_.end()) { + const auto& g_coords = domainInfo_.idxToGridCoord(gridIdx); + auto result = calcPlicOrPlic3CellClass(domainInfo_, g_coords, vofData_, eps_, numIterations_); + it = cache_.insert({gridIdx, std::move(result)}).first; + } + return it->second; + } + + inline const PlicCellClassResult& get(const vec3& pos) { + return get(domainInfo_.posToIdxEx(pos)); + } + + private: + const DomainInfo& domainInfo_; + const VofData& vofData_; + double eps_; + std::size_t numIterations_; + std::unordered_map cache_; + }; + +} // namespace VofFlow diff --git a/src/VofDataUtils/Plic/Plic.cpp b/src/VofDataUtils/Plic/Plic.cpp new file mode 100644 index 0000000..39ee243 --- /dev/null +++ b/src/VofDataUtils/Plic/Plic.cpp @@ -0,0 +1,84 @@ +#include "Plic.h" + +#include +#include + +#include + +#include "../Grid/Gradient.h" +#include "../Misc/Profiling.h" + +VofFlow::PlicCellResult VofFlow::calcPlicCell(const K::Point_3& minP, const K::Point_3& maxP, double f, + const K::Vector_3& gradient, double eps, std::size_t numIterations) { + ZoneScoped; + + // Corners of different phase + const K::Point_3 corner1(CGAL::is_negative(gradient.x()) ? minP.x() : maxP.x(), + CGAL::is_negative(gradient.y()) ? minP.y() : maxP.y(), CGAL::is_negative(gradient.z()) ? minP.z() : maxP.z()); + const K::Point_3 corner0(CGAL::is_negative(gradient.x()) ? maxP.x() : minP.x(), + CGAL::is_negative(gradient.y()) ? maxP.y() : minP.y(), CGAL::is_negative(gradient.z()) ? maxP.z() : minP.z()); + + // Cell + const K::Iso_cuboid_3 cell(minP, maxP, 0); + const double cell_volume = CGAL::to_double(cell.volume()); + + // Normal + const auto squaredLength = gradient.squared_length(); + const K::Vector_3 normal = !CGAL::is_zero(squaredLength) ? -gradient / CGAL::sqrt(squaredLength) : gradient; + + // Diagonal + const K::Vector_3 diag = corner0 - corner1; + + // Iso + double iso_val = f; + double min_val = eps; + double max_val = 1.0 - eps; + PolyCell poly_cell; + std::size_t iter = 0; + double error = std::numeric_limits::max(); + + while (iter < numIterations && error >= eps) { + const K::Point_3 upPoint(corner1 + iso_val * diag); + + poly_cell = PolyCell(cell, K::Plane_3(upPoint, normal)); + + double volume = poly_cell.volume() / cell_volume; + if (volume > f) { + max_val = iso_val; + } else { + min_val = iso_val; + } + + iso_val = (max_val + min_val) * 0.5; + + iter++; + error = std::abs(volume - f); + } + + return {poly_cell, iter, error, corner1, corner0}; +} + +VofFlow::PlicCellResult VofFlow::calcPlicCell(const DomainInfo& domainInfo, const gridCoords_t& g_coords, + const vtkSmartPointer& data, const K::Vector_3& gradient, double eps, std::size_t numIterations) { + ZoneScoped; + + // F + const auto& idx = domainInfo.gridCoordToIdx(g_coords); + const double f = data->GetComponent(idx, 0); + + // PLIC + const auto cellPMin = toPoint_3(domainInfo.coords(g_coords)); + const auto cellPMax = toPoint_3(domainInfo.coords(g_coords[0] + 1, g_coords[1] + 1, g_coords[2] + 1)); + + return calcPlicCell(cellPMin, cellPMax, f, gradient, eps, numIterations); +} + +VofFlow::PlicCellResult VofFlow::calcPlicCell(const DomainInfo& domainInfo, const gridCoords_t& g_coords, + const vtkSmartPointer& data, double eps, std::size_t numIterations) { + ZoneScoped; + + // Gradient + const auto grad = toVector_3(gradient(domainInfo, g_coords, data)); + + return calcPlicCell(domainInfo, g_coords, data, grad, eps, numIterations); +} diff --git a/src/VofDataUtils/Plic/Plic.h b/src/VofDataUtils/Plic/Plic.h new file mode 100644 index 0000000..8a843fd --- /dev/null +++ b/src/VofDataUtils/Plic/Plic.h @@ -0,0 +1,46 @@ +#pragma once + +#include + +#include +#include +#include +#include + +#include "../Grid/DomainInfo.h" +#include "../Grid/GridTypes.h" +#include "../Misc/CgalUtil.h" +#include "PolyCell.h" +#include "Polygon3D.h" + +namespace VofFlow { + struct PlicCellResult { + PolyCell cell; + std::size_t iter = 0; + double err = 0.0; + K::Point_3 baseCorner; + K::Point_3 oppositeCorner; + }; + + struct PlicPolyResult { + Polygon3D poly; + std::size_t iter = 0; + double err = 0.0; + }; + + PlicCellResult calcPlicCell(const K::Point_3& minP, const K::Point_3& maxP, double f, const K::Vector_3& gradient, + double eps = 1e-6, std::size_t numIterations = 20); + + PlicCellResult calcPlicCell(const DomainInfo& domainInfo, const gridCoords_t& g_coords, + const vtkSmartPointer& data, const K::Vector_3& gradient, double eps = 1e-6, + std::size_t numIterations = 20); + + PlicCellResult calcPlicCell(const DomainInfo& domainInfo, const gridCoords_t& g_coords, + const vtkSmartPointer& data, double eps = 1e-6, std::size_t numIterations = 20); + + inline PlicPolyResult calcPlicPoly(const DomainInfo& domainInfo, const gridCoords_t& g_coords, + const vtkSmartPointer& data, double eps = 1e-6, std::size_t numIterations = 20) { + auto tmp = calcPlicCell(domainInfo, g_coords, data, eps, numIterations); + return {tmp.cell.getPolygon(0), tmp.iter, tmp.err}; + } +} // namespace VofFlow diff --git a/src/VofDataUtils/Plic/Plic3.cpp b/src/VofDataUtils/Plic/Plic3.cpp new file mode 100644 index 0000000..25e6e28 --- /dev/null +++ b/src/VofDataUtils/Plic/Plic3.cpp @@ -0,0 +1,184 @@ +#include "Plic3.h" + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "../Misc/Profiling.h" + +namespace { + inline bool hasNonZeroNeighbor(const VofFlow::DomainInfo& domainInfo, const VofFlow::gridCoords_t& g_coords, + const vtkSmartPointer& data, double eps) { + ZoneScoped; + + const auto& dims = domainInfo.cellDims(); + + for (int i_z = std::max(0, g_coords[2] - 1); i_z < std::min(dims[2] - 1, g_coords[2] + 1); i_z++) { + for (int i_y = std::max(0, g_coords[1] - 1); i_y < std::min(dims[1] - 1, g_coords[1] + 1); i_y++) { + for (int i_x = std::max(0, g_coords[0] - 1); i_x < std::min(dims[0] - 1, g_coords[0] + 1); i_x++) { + if (i_x == g_coords[0] && i_y == g_coords[1] && i_z == g_coords[2]) { + continue; + } + const auto& idx = domainInfo.gridCoordToIdx(i_x, i_y, i_z); + if (data->GetComponent(idx, 0) > eps) { + return true; + } + } + } + } + return false; + } +} // namespace + +VofFlow::PlicCellResult VofFlow::calcPlic3Cell(const PolyCell& f3_cell, double f, const K::Vector_3& normal, double eps, + std::size_t numIterations) { + ZoneScoped; + + // poly cell + auto f_cell = f3_cell; + f_cell.invertPlanes(); + + const auto& poly_points = f_cell.polyhedronPoints(); + + // determine corner and opposite corner + const K::Plane_3 distPlane(poly_points[0], normal); + double minSqDist = 0; + double maxSqDist = 0; + std::size_t liquid_corner_idx = 0; + std::size_t vaporous_corner_idx = 0; + + for (std::size_t i = 0; i < poly_points.size(); i++) { + const auto& point = poly_points[i]; + double sqDist = CGAL::to_double(CGAL::squared_distance(distPlane, point)); + if (distPlane.has_on_negative_side(point)) { + sqDist = -sqDist; + } + if (sqDist < minSqDist) { + minSqDist = sqDist; + liquid_corner_idx = i; + } + if (sqDist > maxSqDist) { + maxSqDist = sqDist; + vaporous_corner_idx = i; + } + } + + // Get corner point and normal + const auto corner = poly_points[liquid_corner_idx]; + const auto opposite_corner = poly_points[vaporous_corner_idx]; + + // Iso + const K::Vector_3 poly_cell_diag = opposite_corner - corner; + const auto poly_cell_volume = f_cell.volume(); + const auto total_cell_volume = CGAL::to_double(f_cell.cellCube().volume()); + + double min_val = eps; + double max_val = 1.0 - eps; + double iso_val = std::clamp(f * total_cell_volume / poly_cell_volume, min_val, max_val); + + PolyCell poly_cell; + std::size_t iter = 0; + double error = std::numeric_limits::max(); + + while (iter < numIterations && error >= eps) { + const K::Point_3 upPoint(corner + iso_val * poly_cell_diag); + + poly_cell = PolyCell(f_cell, K::Plane_3(upPoint, normal)); + + double volume = poly_cell.volume() / total_cell_volume; + + if (volume > f) { + max_val = iso_val; + } else { + min_val = iso_val; + } + + iso_val = (max_val + min_val) * 0.5; + + iter++; + error = std::abs(volume - f); + } + + return {poly_cell, iter, error, corner, opposite_corner}; +} + +VofFlow::PlicCellClassResult VofFlow::calcPlic3CellClass(const DomainInfo& domainInfo, const gridCoords_t& g_coords, + const VofData& vofData, double eps, std::size_t numIterations) { + ZoneScoped; + + const double eps_one = 1.0 - eps; + + const vtkIdType idx = domainInfo.gridCoordToIdx(g_coords); + + const double f1 = vofData.vof1st->GetComponent(idx, 0); + const double f2 = vofData.vof2nd->GetComponent(idx, 0); + + // empty cells + if (f1 < eps && f2 < eps) { + return {CellClass::EMPTY, std::nullopt, std::nullopt}; + } + + // full solid cell (interface between full and empty cell is ignored here) + if (f1 > eps_one) { + return {CellClass::FULL_PHASE1, std::nullopt, std::nullopt}; + } + + // full fluid cell (interface between full and empty cell is ignored here) + if (f2 > eps_one) { + return {CellClass::FULL_PHASE2, std::nullopt, std::nullopt}; + } + + // only solid + fluid + if (f1 + f2 > eps_one) { + auto plic = calcPlicCell(domainInfo, g_coords, vofData.vof1st, eps, numIterations); + return {CellClass::INTERFACE_PH1_PH2, plic, std::nullopt}; + } + + // only solid interface + if (f2 < eps) { + auto plic = calcPlicCell(domainInfo, g_coords, vofData.vof1st, eps, numIterations); + return {CellClass::INTERFACE_PHASE1, plic, std::nullopt}; + } + + // only fluid interface + if (f1 < eps) { + if (!hasNonZeroNeighbor(domainInfo, g_coords, vofData.vof1st, eps)) { + // without solid neighbor + auto plic = calcPlicCell(domainInfo, g_coords, vofData.vof2nd, eps, numIterations); + return {CellClass::INTERFACE_PHASE2, plic, std::nullopt}; + } else { + // neighbor cell with solid (Three-Phase - Typ 2) + std::array norm{}; + vofData.vof2ndNorm->GetTuple(idx, norm.data()); + const K::Vector_3 grad{-norm[0], -norm[1], -norm[2]}; + + auto plic = calcPlicCell(domainInfo, g_coords, vofData.vof2nd, grad, eps, numIterations); + return {CellClass::INTERFACE_PHASE2, plic, std::nullopt}; + } + } + + // three-phase cell (Three-Phase - Typ 1 or 3) + if (f1 >= eps && f2 >= eps && f1 + f2 <= eps_one) { + auto f1_result = calcPlicCell(domainInfo, g_coords, vofData.vof1st, eps, numIterations); + + std::array norm{}; + vofData.vof2ndNorm->GetTuple(idx, norm.data()); + const K::Vector_3 normal{norm[0], norm[1], norm[2]}; + + auto f2_result = calcPlic3Cell(f1_result.cell, f2, normal, eps, numIterations); + + return {CellClass::INTERFACE_ALL, f1_result, f2_result}; + } + + // should never reach here! + throw std::runtime_error("Invalid PLIC3!"); +} diff --git a/src/VofDataUtils/Plic/Plic3.h b/src/VofDataUtils/Plic/Plic3.h new file mode 100644 index 0000000..e61ca6a --- /dev/null +++ b/src/VofDataUtils/Plic/Plic3.h @@ -0,0 +1,43 @@ +#pragma once + +#include +#include + +#include "../Grid/DomainInfo.h" +#include "../Grid/GridTypes.h" +#include "../Misc/CgalUtil.h" +#include "../Misc/VofData.h" +#include "Plic.h" +#include "PolyCell.h" + +namespace VofFlow { + enum class CellClass { + EMPTY = 0, + FULL_PHASE1, // Phase 1 = Solid or main phase with regular gradient plic. + FULL_PHASE2, // Phase 2 = Fluid or second phase with polyhedron plic. + INTERFACE_PH1_PH2, + INTERFACE_PHASE1, + INTERFACE_PHASE2, + INTERFACE_ALL, + }; + + struct PlicCellClassResult { + CellClass cellClass = CellClass::EMPTY; + std::optional plic1; + std::optional plic2; + + [[nodiscard]] inline const K::Plane_3& plane1() const { + return plic1->cell.getPlane(0); + } + + [[nodiscard]] inline const K::Plane_3& plane2() const { + return plic2->cell.getPlane(1); + } + }; + + PlicCellResult calcPlic3Cell(const PolyCell& f3_cell, double f, const K::Vector_3& normal, double eps = 1e-6, + std::size_t numIterations = 20); + + PlicCellClassResult calcPlic3CellClass(const DomainInfo& domainInfo, const gridCoords_t& g_coords, + const VofData& vofData, double eps = 1e-6, std::size_t numIterations = 20); +} // namespace VofFlow diff --git a/src/VofDataUtils/Plic/PlicPolyData.cpp b/src/VofDataUtils/Plic/PlicPolyData.cpp new file mode 100644 index 0000000..94286a0 --- /dev/null +++ b/src/VofDataUtils/Plic/PlicPolyData.cpp @@ -0,0 +1,59 @@ +#include "PlicPolyData.h" + +#include +#include + +#include + +#include +#include + +#include "../VtkUtil/VtkUtilArray.h" + +VofFlow::vtkPlicPolyData::vtkPlicPolyData() { + points_ = vtkSmartPointer::New(); + polys_ = vtkSmartPointer::New(); + iterations_ = createVtkArray("Iterations", 1); + error_ = createVtkArray("Error", 1); +} + +void VofFlow::vtkPlicPolyData::addPolygon(const PlicPolyResult& plicResult) { + auto numPointsBefore = points_->GetNumberOfPoints(); + const auto& cellPoints = plicResult.poly.points(); + for (const auto& p : cellPoints) { + points_->InsertNextPoint(CGAL::to_double(p.x()), CGAL::to_double(p.y()), CGAL::to_double(p.z())); + } + std::vector cell(cellPoints.size()); + std::iota(cell.begin(), cell.end(), numPointsBefore); + polys_->InsertNextCell(static_cast(cell.size()), cell.data()); + iterations_->InsertNextValue(static_cast(plicResult.iter)); + error_->InsertNextValue(static_cast(plicResult.err)); +} + +vtkSmartPointer VofFlow::vtkPlicPolyData::getPolyData() const { + vtkSmartPointer polyData = vtkSmartPointer::New(); + polyData->SetPoints(points_); + polyData->SetPolys(polys_); + polyData->GetCellData()->AddArray(iterations_); + polyData->GetCellData()->AddArray(error_); + return polyData; +} + +VofFlow::vtkPlic3PolyData::vtkPlic3PolyData() { + cellClass_ = createVtkArray("CellClass", 1); + phaseType_ = createVtkArray("PhaseType", 1); +} + +void VofFlow::vtkPlic3PolyData::addPolygon(VofFlow::CellClass cellClass, unsigned char phaseType, + const VofFlow::PlicPolyResult& plicResult) { + vtkPlic_.addPolygon(plicResult); + cellClass_->InsertNextValue(static_cast(cellClass)); + phaseType_->InsertNextValue(static_cast(phaseType)); +} + +vtkSmartPointer VofFlow::vtkPlic3PolyData::getPolyData() const { + auto polyData = vtkPlic_.getPolyData(); + polyData->GetCellData()->AddArray(cellClass_); + polyData->GetCellData()->AddArray(phaseType_); + return polyData; +} diff --git a/src/VofDataUtils/Plic/PlicPolyData.h b/src/VofDataUtils/Plic/PlicPolyData.h new file mode 100644 index 0000000..21d916c --- /dev/null +++ b/src/VofDataUtils/Plic/PlicPolyData.h @@ -0,0 +1,43 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include + +#include "Plic.h" +#include "Plic3.h" + +namespace VofFlow { + class vtkPlicPolyData { + public: + vtkPlicPolyData(); + + void addPolygon(const PlicPolyResult& plicResult); + + [[nodiscard]] vtkSmartPointer getPolyData() const; + + protected: + vtkSmartPointer points_; + vtkSmartPointer polys_; + vtkSmartPointer iterations_; + vtkSmartPointer error_; + }; + + class vtkPlic3PolyData { + public: + vtkPlic3PolyData(); + + void addPolygon(CellClass cellClass, unsigned char phaseType, const PlicPolyResult& plicResult); + + [[nodiscard]] vtkSmartPointer getPolyData() const; + + protected: + vtkPlicPolyData vtkPlic_; + vtkSmartPointer cellClass_; + vtkSmartPointer phaseType_; + }; +} // namespace VofFlow diff --git a/src/VofDataUtils/Plic/PlicUtil.cpp b/src/VofDataUtils/Plic/PlicUtil.cpp new file mode 100644 index 0000000..4993e47 --- /dev/null +++ b/src/VofDataUtils/Plic/PlicUtil.cpp @@ -0,0 +1,52 @@ +#include "PlicUtil.h" + +#include "../Misc/Profiling.h" + +VofFlow::PlicCellClassResult VofFlow::calcPlicOrPlic3CellClass(const DomainInfo& domainInfo, + const gridCoords_t& g_coords, const VofData& vofData, double eps, std::size_t numIterations) { + ZoneScoped; + + if (vofData.vof2nd != nullptr) { + return calcPlic3CellClass(domainInfo, g_coords, vofData, eps, numIterations); + } else { + const auto& idx = domainInfo.gridCoordToIdx(g_coords); + const double f = vofData.vof1st->GetComponent(idx, 0); + if (f < eps) { + return {CellClass::EMPTY, std::nullopt, std::nullopt}; + } else if (f > 1.0 - eps) { + return {CellClass::FULL_PHASE1, std::nullopt, std::nullopt}; + } else { + auto plic = calcPlicCell(domainInfo, g_coords, vofData.vof1st, eps); + return {CellClass::INTERFACE_PHASE1, plic, std::nullopt}; + } + } +} + +unsigned char VofFlow::getPhaseIdx(const PlicCellClassResult& plicResult, const vec3& pos) { + ZoneScoped; + + switch (plicResult.cellClass) { + case CellClass::EMPTY: + return 0; + case CellClass::FULL_PHASE1: + return 1; + case CellClass::FULL_PHASE2: + return 2; + case CellClass::INTERFACE_PH1_PH2: + return !plicResult.plane1().has_on_positive_side(toPoint_3(pos)) ? 1 : 2; + case CellClass::INTERFACE_PHASE1: + return plicResult.plane1().has_on_positive_side(toPoint_3(pos)) ? 0 : 1; + case CellClass::INTERFACE_PHASE2: + return plicResult.plane1().has_on_positive_side(toPoint_3(pos)) ? 0 : 2; + case CellClass::INTERFACE_ALL: + const auto p = toPoint_3(pos); + if (!plicResult.plane1().has_on_positive_side(p)) { + return 1; + } else if (!plicResult.plane2().has_on_positive_side(p)) { + return 2; + } else { + return 0; + } + } + return 0; // Should not happen +} diff --git a/src/VofDataUtils/Plic/PlicUtil.h b/src/VofDataUtils/Plic/PlicUtil.h new file mode 100644 index 0000000..ce966c4 --- /dev/null +++ b/src/VofDataUtils/Plic/PlicUtil.h @@ -0,0 +1,17 @@ +#pragma once + +#include + +#include "../Grid/DomainInfo.h" +#include "../Grid/GridTypes.h" +#include "../Math/Vector.h" +#include "../Misc/VofData.h" +#include "Plic.h" +#include "Plic3.h" + +namespace VofFlow { + PlicCellClassResult calcPlicOrPlic3CellClass(const DomainInfo& domainInfo, const gridCoords_t& g_coords, + const VofData& vofData, double eps = 1e-6, std::size_t numIterations = 20); + + unsigned char getPhaseIdx(const PlicCellClassResult& plicResult, const vec3& pos); +} // namespace VofFlow diff --git a/src/VofDataUtils/Plic/PolyCell.cpp b/src/VofDataUtils/Plic/PolyCell.cpp new file mode 100644 index 0000000..2fb1d52 --- /dev/null +++ b/src/VofDataUtils/Plic/PolyCell.cpp @@ -0,0 +1,178 @@ +#include "PolyCell.h" + +#include +#include +#include + +#include +#include +#include + +#include "../Misc/Profiling.h" + +namespace { + using VofFlow::K; + + struct Intersection_visitor : public boost::static_visitor> { + std::vector operator()(const K::Point_3& p) const { + return {p}; + } + std::vector operator()(const K::Segment_3& s) const { + return {s.source(), s.target()}; + } + std::vector operator()(const K::Triangle_3& t) const { + return {t.vertex(0), t.vertex(1), t.vertex(2)}; + } + std::vector operator()(const std::vector& v) const { + return v; + } + }; + + inline double polyhedronVolume(const std::vector& points) { + double volume = 0.0; + CGAL::Delaunay_triangulation_3 delaunay(points.begin(), points.end()); + for (const auto& cell : delaunay.finite_cell_handles()) { + volume += CGAL::to_double(delaunay.tetrahedron(cell).volume()); + } + return volume; + } +} // namespace + +VofFlow::PolyCell::PolyCell() : cell_() {} + +VofFlow::PolyCell::PolyCell(const K::Iso_cuboid_3& cell) : cell_(cell) { + ZoneScoped; + + updatePolyhedronPoints(); +} + +VofFlow::PolyCell::PolyCell(const K::Iso_cuboid_3& cell, const K::Plane_3& intersect_plane) : cell_(cell) { + ZoneScoped; + + auto p = intersectPlane(intersect_plane); + planes_.emplace_back(std::move(p)); + + updatePolyhedronPoints(); +} + +VofFlow::PolyCell::PolyCell(const PolyCell& cell, const K::Plane_3& intersect_plane) + : cell_(cell.cell_), + planes_(cell.planes_) { + ZoneScoped; + + auto p = intersectPlane(intersect_plane); + planes_.emplace_back(std::move(p)); + + updatePolyhedronPoints(); +} + +VofFlow::Polygon3D VofFlow::PolyCell::getPolygon(std::size_t i) const { + ZoneScoped; + + const auto& p = planes_.at(i); + return {p.plane, p.intersections}; +} + +double VofFlow::PolyCell::volume() const { + ZoneScoped; + + return polyhedronVolume(polyhedron_points_); +} + +void VofFlow::PolyCell::invertPlanes() { + ZoneScoped; + + for (auto& p : planes_) { + p.plane = p.plane.opposite(); + } + updatePolyhedronPoints(); +} + +VofFlow::PolyCell::Plane VofFlow::PolyCell::intersectPlane(const K::Plane_3& intersect_plane) { + ZoneScoped; + + Plane p{intersect_plane, {}}; + + const auto result = intersection(cell_, p.plane); + if (result) { + p.intersections = boost::apply_visitor(Intersection_visitor(), *result); + } + return p; +} + +void VofFlow::PolyCell::updatePolyhedronPoints() { + ZoneScoped; + + polyhedron_points_.clear(); + if (planes_.empty()) { + polyhedron_points_.reserve(8); + for (int v = 0; v < 8; v++) { + polyhedron_points_.push_back(cell_.vertex(v)); + } + } else if (planes_.size() == 1) { + polyhedron_points_ = planes_[0].intersections; + for (int v = 0; v < 8; v++) { + const auto& vert = cell_.vertex(v); + if (planes_[0].plane.has_on_negative_side(vert)) { + polyhedron_points_.push_back(vert); + } + } + } else { + auto addPointIfValid = [this](const K::Point_3& point, const Plane* const ignore = nullptr) { + if (std::all_of(planes_.cbegin(), planes_.cend(), [&point, &ignore](const auto& plane) { + return &plane == ignore || plane.plane.has_on_negative_side(point); + })) { + polyhedron_points_.push_back(point); + } + }; + + // Add intersection points from all planes + for (const auto& plane : planes_) { + for (const auto& point : plane.intersections) { + addPointIfValid(point, &plane); + } + } + // Add cube vertices + for (int v = 0; v < 8; v++) { + addPointIfValid(cell_.vertex(v)); + } + + // Add intersections between planes + // TODO only implemented for 2 planes + if (planes_.size() > 2) { + throw std::runtime_error("More than 2 planes not implemented!"); + } + + const auto result = CGAL::intersection(planes_[0].plane, planes_[1].plane); + if (result) { + if (const K::Line_3* l = boost::get(&*result)) { + const auto result2 = CGAL::intersection(cell_, *l); + if (result2) { + if (const K::Segment_3* s = boost::get(&*result2)) { + polyhedron_points_.push_back(s->source()); + polyhedron_points_.push_back(s->target()); + } else { + const K::Point_3* p = boost::get(&*result2); + polyhedron_points_.push_back(*p); + } + } + } else { + // Planes are identical. + // `addPointIfValid` checks for not on plane, no points are added so far. Just add points from plane 1. + polyhedron_points_.insert(polyhedron_points_.end(), planes_[0].intersections.begin(), + planes_[0].intersections.end()); + } + } else { + // Parallel planes + // Case 1 - normals pointing away from each other: + // All points are passing the `addPointIfValid` test. Nothing to do here. + // Case 2 - normals pointing in the same direction: + // The foremost plane is outside the volumed defined by the other plane. The points of the other plane are + // passing `addPointIfValid`, the points of the foremost plane not. Nothing to do here. + // Case 3 - normals are pointing towards each other: + // This means the second plane will define the volume spanned by the first plane completely as outside. The + // empty point set is the correct solution. Both planes remove the points from the other plane in + // `addPointIfValid`. Nothing to do here. + } + } +} diff --git a/src/VofDataUtils/Plic/PolyCell.h b/src/VofDataUtils/Plic/PolyCell.h new file mode 100644 index 0000000..8e665a8 --- /dev/null +++ b/src/VofDataUtils/Plic/PolyCell.h @@ -0,0 +1,56 @@ +#pragma once + +#include +#include + +#include +#include +#include + +#include "../Misc/CgalUtil.h" +#include "Polygon3D.h" + +namespace VofFlow { + /* + * A PolyCell is defined by an iso cuboid and intersection planes. + */ + class PolyCell { + public: + PolyCell(); + explicit PolyCell(const K::Iso_cuboid_3& cell); + PolyCell(const K::Iso_cuboid_3& cell, const K::Plane_3& intersect_plane); + PolyCell(const PolyCell& cell, const K::Plane_3& intersect_plane); + + [[nodiscard]] const K::Iso_cuboid_3& cellCube() const { + return cell_; + } + + [[nodiscard]] const std::vector& polyhedronPoints() const { + return polyhedron_points_; + } + + [[nodiscard]] inline const K::Plane_3& getPlane(std::size_t i) const { + return planes_.at(i).plane; + } + + [[nodiscard]] Polygon3D getPolygon(std::size_t i) const; + + [[nodiscard]] double volume() const; + + void invertPlanes(); + + protected: + struct Plane { + K::Plane_3 plane; + std::vector intersections; // only the intersections with the cube + }; + + Plane intersectPlane(const K::Plane_3& intersect_plane); + + void updatePolyhedronPoints(); + + K::Iso_cuboid_3 cell_; + std::vector planes_; + std::vector polyhedron_points_; + }; +} // namespace VofFlow diff --git a/src/VofDataUtils/Plic/Polygon3D.cpp b/src/VofDataUtils/Plic/Polygon3D.cpp new file mode 100644 index 0000000..375f8a4 --- /dev/null +++ b/src/VofDataUtils/Plic/Polygon3D.cpp @@ -0,0 +1,90 @@ +#include "Polygon3D.h" + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "../Misc/Profiling.h" + +VofFlow::Polygon3D::Polygon3D(const K::Plane_3& plane, const std::vector& points) : plane_(plane) { + ZoneScoped; + + // Project points on plain + // Note: to_2d() result coordinate system seems to be based on the planes normal length. If normal is normalized, + // coords seems to be at the same scale as the 3D space. + std::vector points_2d(points.size()); + for (std::size_t i = 0; i < points.size(); i++) { + points_2d[i] = plane_.to_2d(points[i]); + } + + // Calc Convex hull to make sure points are unique (in case of PLIC3 vertices from multiple tetrahedron could be + // doubled), also for sorting. + CGAL::convex_hull_2(points_2d.begin(), points_2d.end(), std::back_inserter(points_2d_)); + + if (points_2d_.size() < 3) { + throw std::runtime_error("Polygon validation failed!"); + } + + // Reproject ordered points back to 3d + points_3d_.resize(points_2d_.size()); + for (std::size_t i = 0; i < points_2d_.size(); i++) { + points_3d_[i] = plane_.to_3d(points_2d_[i]); + } +} + +bool VofFlow::Polygon3D::isIn2DPolygon(const K::Point_2& p) const { + ZoneScoped; + + return CGAL::bounded_side_2(points_2d_.begin(), points_2d_.end(), p) != CGAL::ON_UNBOUNDED_SIDE; +} + +bool VofFlow::Polygon3D::isOrthogonalProjectionInsidePolygon(const K::Point_3& p) const { + ZoneScoped; + + const auto pOnPlane = plane_.to_2d(p); + return isIn2DPolygon(pOnPlane); +} + +VofFlow::K::Point_3 VofFlow::Polygon3D::findNearestPoint(const K::Point_3& p) const { + ZoneScoped; + + const auto pOnPlane = plane_.to_2d(p); + if (isIn2DPolygon(pOnPlane)) { + return plane_.to_3d(pOnPlane); + } + // Iterate over all neighboring pairs of points and find the closest point on the corresponding line segment. + double minDist = std::numeric_limits::max(); + K::Point_2 bestP; + for (std::size_t i = 0; i < points_2d_.size(); i++) { + const auto& p1 = points_2d_[i]; + const auto& p2 = points_2d_[(i + 1) % points_2d_.size()]; + K::Segment_2 segment(p1, p2); + const K::Point_2 pOnLine = segment.supporting_line().projection(pOnPlane); + + const auto dist_toP1 = CGAL::squared_distance(p1, pOnLine); + const auto dist_toP2 = CGAL::squared_distance(p2, pOnLine); + const auto dist_P1P2 = CGAL::squared_distance(p1, p2); + K::Point_2 nearestP; + if (dist_toP1 <= dist_P1P2 && dist_toP2 <= dist_P1P2) { // point in between + nearestP = pOnLine; + } else if (dist_toP1 < dist_toP2) { // not in between, nearer to P1 + nearestP = p1; + } else { // not in between, nearer to P2 + nearestP = p2; + } + + const auto dist = CGAL::squared_distance(nearestP, pOnPlane); + if (dist < minDist) { + minDist = CGAL::to_double(dist); + bestP = nearestP; + } + } + + return plane_.to_3d(bestP); +} diff --git a/src/VofDataUtils/Plic/Polygon3D.h b/src/VofDataUtils/Plic/Polygon3D.h new file mode 100644 index 0000000..f921520 --- /dev/null +++ b/src/VofDataUtils/Plic/Polygon3D.h @@ -0,0 +1,44 @@ +#pragma once + +#include + +#include +#include +#include +#include + +#include "../Misc/CgalUtil.h" + +namespace VofFlow { + class Polygon3D { + public: + Polygon3D(const K::Plane_3& plane, const std::vector& points); + + [[nodiscard]] const K::Plane_3& plane() const { + return plane_; + } + + [[nodiscard]] const std::vector& points() const { + return points_3d_; + } + + [[nodiscard]] const std::vector& points2() const { + return points_2d_; + } + + [[nodiscard]] CGAL::Bbox_2 bbox2() const { + return CGAL::bbox_2(points_2d_.begin(), points_2d_.end()); + } + + [[nodiscard]] bool isIn2DPolygon(const K::Point_2& p) const; + + [[nodiscard]] bool isOrthogonalProjectionInsidePolygon(const K::Point_3& p) const; + + [[nodiscard]] K::Point_3 findNearestPoint(const K::Point_3& p) const; + + protected: + K::Plane_3 plane_; + std::vector points_3d_; + std::vector points_2d_; + }; +} // namespace VofFlow diff --git a/src/VofDataUtils/VtkUtil/VtkUtilArray.cpp b/src/VofDataUtils/VtkUtil/VtkUtilArray.cpp new file mode 100644 index 0000000..f339b4a --- /dev/null +++ b/src/VofDataUtils/VtkUtil/VtkUtilArray.cpp @@ -0,0 +1,25 @@ +#include "VtkUtilArray.h" + +#include +#include + +namespace { + struct ShallowCopyWorker { + vtkSmartPointer copy_; + + template + void operator()(vtkAOSDataArrayTemplate* dataArray) { + copy_ = vtkSmartPointer>::New(); + copy_->ShallowCopy(dataArray); + } + }; +} // namespace + +vtkSmartPointer VofFlow::shallowCopy(const vtkSmartPointer& data) { + ShallowCopyWorker worker; + typedef vtkArrayDispatch::DispatchByValueType Dispatcher; + if (!Dispatcher::Execute(data, worker)) { + throw std::runtime_error("Cannot dispatch array worker!"); + } + return worker.copy_; +} diff --git a/src/VofDataUtils/VtkUtil/VtkUtilArray.h b/src/VofDataUtils/VtkUtil/VtkUtilArray.h new file mode 100644 index 0000000..a095580 --- /dev/null +++ b/src/VofDataUtils/VtkUtil/VtkUtilArray.h @@ -0,0 +1,104 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "../Math/Vector.h" + +namespace VofFlow { + vtkSmartPointer shallowCopy(const vtkSmartPointer& data); + + template + inline vtkSmartPointer createVtkArray(const char* name, int numComponents = 1, vtkIdType numTuples = 0) { + vtkSmartPointer a = vtkSmartPointer::New(); + a->SetName(name); + a->SetNumberOfComponents(numComponents); + if (numTuples > 0) { + a->SetNumberOfTuples(numTuples); + } + return a; + } + + template + inline vtkSmartPointer createVtkArray(const char* name, int numComponents, vtkIdType numTuples, + double fill) { + auto a = createVtkArray(name, numComponents, numTuples); + a->Fill(fill); + return a; + } + + template + inline vtkSmartPointer createVtkArray(const char* name, + const std::vector& data) { + auto a = createVtkArray(name, 1, static_cast(data.size())); + std::memcpy(a->GetVoidPointer(0), data.data(), sizeof(typename ArrayT::ValueType) * data.size()); + return a; + } + + inline vtkSmartPointer createVtkArray(const char* name, const std::vector& data) { + auto a = createVtkArray(name, 3, static_cast(data.size())); + std::memcpy(a->GetVoidPointer(0), data.data(), sizeof(vec3) * data.size()); + return a; + } + + inline vtkSmartPointer createVtkStringArrayFromString(const char* name, const std::string& str) { + vtkSmartPointer a = vtkSmartPointer::New(); + a->SetName(name); + a->SetNumberOfValues(1); + a->SetValue(0, str); + return a; + } + + template + void appendArray(ArrayT* a, ArrayT* b) { + if (a->GetNumberOfComponents() != b->GetNumberOfComponents()) { + throw std::runtime_error("Bad array format!"); + } + auto* dest = a->WritePointer(a->GetNumberOfValues(), b->GetNumberOfValues()); + const auto* src = b->GetPointer(0); + std::memcpy(dest, src, b->GetNumberOfValues() * sizeof(typename ArrayT::ValueType)); + } + + template + void appendArray(vtkSmartPointer& a, vtkSmartPointer& b) { + appendArray(a.Get(), b.Get()); + } + + template + inline typename ArrayT::ValueType* getTuplePointer(ArrayT* array, vtkIdType tupleIdx) { + return array->GetPointer(tupleIdx * array->GetNumberOfComponents()); + } + + template + inline typename ArrayT::ValueType* getTuplePointer(const vtkSmartPointer& array, vtkIdType tupleIdx) { + return getTuplePointer(array.Get(), tupleIdx); + } + + template + std::vector getArrayValues(vtkDataArray* data) { + const vtkIdType size = data->GetNumberOfTuples(); + std::vector v(size); + for (vtkIdType i = 0; i < size; i++) { + v[i] = static_cast(data->GetComponent(i, 0)); + } + return v; + } + + inline void appendArrayName(vtkDataArray* data, const std::string& nameExtra) { + data->SetName((std::string(data->GetName()) + nameExtra).c_str()); + } + + inline vec3& getVec3TuplePointer(vtkFloatArray* array, vtkIdType tupleIdx) { + // assume NumberOfComponents == 3 + return *reinterpret_cast(getTuplePointer(array, tupleIdx)); + } + +} // namespace VofFlow diff --git a/src/VofDataUtils/VtkUtil/VtkUtilMPI.h b/src/VofDataUtils/VtkUtil/VtkUtilMPI.h new file mode 100644 index 0000000..d0ef73c --- /dev/null +++ b/src/VofDataUtils/VtkUtil/VtkUtilMPI.h @@ -0,0 +1,64 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include + +namespace VofFlow { + template + int AllGatherVSameLength(vtkMPIController* mpiController, const T* sendBuffer, T* recvBuffer, vtkIdType length) { + const int numProcesses = mpiController->GetNumberOfProcesses(); + std::vector recvLengths(numProcesses, length); + std::vector recvOffsets(numProcesses); + for (int i = 0; i < numProcesses; ++i) { + recvOffsets[i] = i * length; + } + return mpiController->AllGatherV(sendBuffer, recvBuffer, length, recvLengths.data(), recvOffsets.data()); + } + + template + inline int AllGatherVSameLength(vtkMPIController* mpiController, const std::vector& sendBuffer, + std::vector& recvBuffer) { + const int numProcesses = mpiController->GetNumberOfProcesses(); + if (sendBuffer.size() * numProcesses != recvBuffer.size()) { + throw std::runtime_error("Bad buffer sizes!"); + } + return AllGatherVSameLength(mpiController, sendBuffer.data(), recvBuffer.data(), sendBuffer.size()); + } + + template + int SendVector(vtkMPIController* mpiController, const std::vector& data, int remoteProcessId, int tag) { + static_assert(std::is_trivially_copyable_v, "Vector type must be trivially copyable!"); + std::size_t size = data.size(); + mpiController->Send(&size, 1, remoteProcessId, tag); + if (size > 0) { + mpiController->Send(reinterpret_cast(data.data()), size * sizeof(T), remoteProcessId, tag); + } + return 1; + } + + template + int ReceiveVector(vtkMPIController* mpiController, std::vector& data, int remoteProcessId, int tag) { + static_assert(std::is_trivially_copyable_v, "Vector type must be trivially copyable!"); + std::size_t size = 0; + mpiController->Receive(&size, 1, remoteProcessId, tag); + if (size > 0) { + data.resize(size); + mpiController->Receive(reinterpret_cast(data.data()), size * sizeof(T), remoteProcessId, tag); + } else { + data.clear(); + } + return 1; + } + + template + T ReduceSum(vtkMPIController* mpiController, T value, int destProcessId = 0) { + T globalValue{}; + mpiController->Reduce(&value, &globalValue, 1, vtkCommunicator::StandardOperations::SUM_OP, destProcessId); + return globalValue; + } +} // namespace VofFlow diff --git a/src/VofDataUtils/VtkUtil/VtkUtilWriter.cpp b/src/VofDataUtils/VtkUtil/VtkUtilWriter.cpp new file mode 100644 index 0000000..b46d649 --- /dev/null +++ b/src/VofDataUtils/VtkUtil/VtkUtilWriter.cpp @@ -0,0 +1,54 @@ +#include "VtkUtilWriter.h" + +#include +#include +#include +#include + +void VofFlow::writeData(vtkPolyData* data, const std::string& path, bool compressed) { + vtkSmartPointer writer = vtkSmartPointer::New(); + writer->SetInputData(data); + writer->SetFileName(path.c_str()); + writer->SetDataModeToAppended(); + writer->SetEncodeAppendedData(false); + if (compressed) { + writer->SetCompressorTypeToZLib(); + writer->SetCompressionLevel(1); + } else { + writer->SetCompressorTypeToNone(); + } + writer->SetHeaderTypeToUInt64(); + writer->Write(); +} + +void VofFlow::writeData(vtkRectilinearGrid* data, const std::string& path, bool compressed) { + vtkSmartPointer writer = vtkSmartPointer::New(); + writer->SetInputData(data); + writer->SetFileName(path.c_str()); + writer->SetDataModeToAppended(); + writer->SetEncodeAppendedData(false); + if (compressed) { + writer->SetCompressorTypeToZLib(); + writer->SetCompressionLevel(1); + } else { + writer->SetCompressorTypeToNone(); + } + writer->SetHeaderTypeToUInt64(); + writer->Write(); +} + +void VofFlow::writeData(vtkImageData* data, const std::string& path, bool compressed) { + vtkSmartPointer writer = vtkSmartPointer::New(); + writer->SetInputData(data); + writer->SetFileName(path.c_str()); + writer->SetDataModeToAppended(); + writer->SetEncodeAppendedData(false); + if (compressed) { + writer->SetCompressorTypeToZLib(); + writer->SetCompressionLevel(1); + } else { + writer->SetCompressorTypeToNone(); + } + writer->SetHeaderTypeToUInt64(); + writer->Write(); +} diff --git a/src/VofDataUtils/VtkUtil/VtkUtilWriter.h b/src/VofDataUtils/VtkUtil/VtkUtilWriter.h new file mode 100644 index 0000000..0f2e13c --- /dev/null +++ b/src/VofDataUtils/VtkUtil/VtkUtilWriter.h @@ -0,0 +1,13 @@ +#pragma once + +#include + +#include +#include +#include + +namespace VofFlow { + void writeData(vtkPolyData* data, const std::string& path, bool compressed = false); + void writeData(vtkRectilinearGrid* data, const std::string& path, bool compressed = false); + void writeData(vtkImageData* data, const std::string& path, bool compressed = false); +} // namespace VofFlow diff --git a/src/VofDataUtils/vtk.module b/src/VofDataUtils/vtk.module new file mode 100644 index 0000000..b3a47b0 --- /dev/null +++ b/src/VofDataUtils/vtk.module @@ -0,0 +1,6 @@ +NAME + VofFlow::VofDataUtils +DEPENDS + VTK::CommonCore + VTK::FiltersCore + VTK::ParallelMPI diff --git a/src/VofTracking/Advection.cpp b/src/VofTracking/Advection.cpp new file mode 100644 index 0000000..aef7fe6 --- /dev/null +++ b/src/VofTracking/Advection.cpp @@ -0,0 +1,428 @@ +#include "Advection.h" + +#include +#include +#include +#include +#include + +#include +#include + +#include "Grid/DataInterpolation.h" +#include "Grid/DomainInfo.h" +#include "Grid/GridIterator.h" +#include "Misc/CgalUtil.h" +#include "Misc/Profiling.h" +#include "Plic/CachedPlic.h" +#include "Plic/PlicUtil.h" + +namespace { + // Returns mirrorInfo_t of which bounds where mirrored. + // Return value should be used in mirrorVelocity after values are read from the mirrored position. + inline VofFlow::mirrorInfo_t mirrorPos(VofFlow::vec3& pos, const VofFlow::DomainInfo& domainInfo, + const VofFlow::mirrorInfo_t& mirrorInfo) { + // Fast return if mirroring is off + if (std::none_of(mirrorInfo.begin(), mirrorInfo.end(), [](bool b) { return b; })) { + return mirrorInfo; + } + + VofFlow::mirrorInfo_t outInfo{}; + const auto& b = domainInfo.globalBounds(); + if (mirrorInfo[0] && pos.x < b[0]) { + pos.x = static_cast(b[0] + (b[0] - pos.x)); + outInfo[0] = true; + } + if (mirrorInfo[1] && pos.x > b[1]) { + pos.x = static_cast(b[1] - (pos.x - b[1])); + outInfo[1] = true; + } + if (mirrorInfo[2] && pos.y < b[2]) { + pos.y = static_cast(b[2] + (b[2] - pos.y)); + outInfo[2] = true; + } + if (mirrorInfo[3] && pos.y > b[3]) { + pos.y = static_cast(b[3] - (pos.y - b[3])); + outInfo[3] = true; + } + if (mirrorInfo[4] && pos.z < b[4]) { + pos.z = static_cast(b[4] + (b[4] - pos.z)); + outInfo[4] = true; + } + if (mirrorInfo[5] && pos.z > b[5]) { + pos.z = static_cast(b[5] - (pos.z - b[5])); + outInfo[5] = true; + } + return outInfo; + } + + inline void mirrorVelocity(VofFlow::vec3& velocity, const VofFlow::mirrorInfo_t& mirrorInfo) { + if (mirrorInfo[0] || mirrorInfo[1]) { + velocity.x = -velocity.x; + } + if (mirrorInfo[2] || mirrorInfo[3]) { + velocity.y = -velocity.y; + } + if (mirrorInfo[4] || mirrorInfo[5]) { + velocity.z = -velocity.z; + } + } + + inline VofFlow::vec3 euler(const VofFlow::vec3& pos0, const VofFlow::DomainInfo& domainInfo, + vtkDataArray* velocityArray0, vtkDataArray* velocityArray1, const double deltaT, const int numSteps, + const VofFlow::mirrorInfo_t& mirrorInfo) { + VofFlow::vec3 pos = pos0; + + const float dt = static_cast(deltaT / numSteps); + const float invNumSteps = 1.0f / static_cast(numSteps); + + for (int i = 0; i < numSteps; ++i) { + // Time in range [0, 1] for linear interpolation between time step data. + const float h = static_cast(i) * invNumSteps; + + const auto vel = interpolateMultiCellDataVec3(velocityArray0, velocityArray1, domainInfo, pos, h); + pos += dt * vel; + + // We only need to mirror the final position, as only the velocity from the original position is used. + mirrorPos(pos, domainInfo, mirrorInfo); + } + + return pos; + } + + inline VofFlow::vec3 rungeKutta4(const VofFlow::vec3& pos0, const VofFlow::DomainInfo& domainInfo, + vtkDataArray* velocityArray0, vtkDataArray* velocityArray1, const double deltaT, const int numSteps, + const VofFlow::mirrorInfo_t& mirrorInfo) { + VofFlow::vec3 pos = pos0; + + const float inv6 = 1.0f / 6.0f; + const float dt = static_cast(deltaT / numSteps); + const float invNumSteps = 1.0f / static_cast(numSteps); + + for (int i = 0; i < numSteps; ++i) { + // Time in range [0, 1] for linear interpolation between time step data. + const float h0 = static_cast(i) * invNumSteps; + const float h05 = (static_cast(i) + 0.5f) * invNumSteps; + const float h1 = static_cast(i + 1) * invNumSteps; + + // k1 + auto k1 = interpolateMultiCellDataVec3(velocityArray0, velocityArray1, domainInfo, pos, h0); + // No mirroring needed, original pos is always in domain. + + // k2 + VofFlow::vec3 k1p = pos + k1 * dt * 0.5f; + const auto mirror_k1 = mirrorPos(k1p, domainInfo, mirrorInfo); + auto k2 = interpolateMultiCellDataVec3(velocityArray0, velocityArray1, domainInfo, k1p, h05); + mirrorVelocity(k2, mirror_k1); + + // k3 + VofFlow::vec3 k2p = pos + k2 * dt * 0.5f; + const auto mirror_k2 = mirrorPos(k2p, domainInfo, mirrorInfo); + auto k3 = interpolateMultiCellDataVec3(velocityArray0, velocityArray1, domainInfo, k2p, h05); + mirrorVelocity(k3, mirror_k2); + + // k4 + VofFlow::vec3 k3p = pos + k3 * dt; + const auto mirror_k3 = mirrorPos(k3p, domainInfo, mirrorInfo); + auto k4 = interpolateMultiCellDataVec3(velocityArray0, velocityArray1, domainInfo, k3p, h1); + mirrorVelocity(k4, mirror_k3); + + // sum + pos += dt * inv6 * (k1 + 2.0f * k2 + 2.0f * k3 + k4); + mirrorPos(pos, domainInfo, mirrorInfo); + } + + return pos; + } + + inline void neighborCorrector(const VofFlow::DomainInfo& domainInfo, + const std::unordered_map>& oldCellToParticleIdMap, + VofFlow::UpdateParticle& updateParticle, const std::vector& particlesPos, + const std::vector& oldParticlesPos, vtkDataArray* vof, double eps) { + ZoneScoped; + + auto oldGridCoord = domainInfo.posToGridCoordEx(oldParticlesPos[updateParticle.listIdx]); + + float minDist = std::numeric_limits::max(); + int bestNeighborIdx = -1; + + // Loop over neighbor cells + for (const VofFlow::gridCoords_t& coords : VofFlow::GridRange(domainInfo.cellDims(), oldGridCoord, 1)) { + // Get particles in neighbor cells + auto it = oldCellToParticleIdMap.find(domainInfo.gridCoordToIdx(coords)); + if (it != oldCellToParticleIdMap.end()) { + // Loop over particles in each cell + for (const auto& neighborIdx : it->second) { + // Ignore the current particle itself + if (updateParticle.listIdx == neighborIdx) { + continue; + } + // Get f value (after advection) of neighbor particle (if in domain) + const auto gridIdx = domainInfo.posToIdx(particlesPos[neighborIdx]); + if (gridIdx.has_value()) { + const double f = vof->GetComponent(gridIdx.value(), 0); + // If cell is filled check which neighbor is closest + if (f > eps) { + const float dist = + lengthSq(oldParticlesPos[neighborIdx] - oldParticlesPos[updateParticle.listIdx]); + if (minDist > dist) { + minDist = dist; + bestNeighborIdx = static_cast(neighborIdx); + } + } + } + } + } + } + if (bestNeighborIdx > -1) { + const auto oldPos = oldParticlesPos[updateParticle.listIdx]; + const auto neighborDisplacement = particlesPos[bestNeighborIdx] - oldParticlesPos[bestNeighborIdx]; + const auto newPos = oldPos + neighborDisplacement; + const auto idx = domainInfo.posToIdx(newPos); + if (idx.has_value()) { + updateParticle.pos = newPos; + updateParticle.gridIdx = idx.value(); + updateParticle.f = vof->GetComponent(idx.value(), 0); + } else { + // Displacement vector of closest neighbor particle applied to current particle goes out of domain + // bounds. Current logic is, do nothing here. This means, if the closest neighbor displacement goes out + // of bounds we do not apply neighbor correction. TODO: Should we move this check to the neighbor search + // above, to find closest neighbor which displacement stays in bounds instead? + // std::cout << "Smart corrector displacement out of bounds!" << std::endl; + } + } + } + + inline VofFlow::vec3 clampCell(const VofFlow::vec3& pos, const VofFlow::vec3& cellMin, + const VofFlow::vec3& cellMax) { + ZoneScoped; + + // VTK ComputeStructuredCoordinates assumes inside as "lower_bound <= p < upper_bound". Therefore, + // clamp to [lower_bound, upper_bound - eps] and use the smallest possible float step as eps. + // In addition, there could be floating point issues from rounding double to float, i.e., if cellMin is + // initially defined as double, the closest float might be slightly smaller. If then + // computeStructuredCoordinates assumes the double value as exact cell boundary, this could technically move + // the position to the wrong cell. To compensate this, an additional epsilon is added to the clamp boundaries: + // [lower_bound + eps, upper_bound - 2 * eps]. + return VofFlow::clamp(pos, VofFlow::nextafter(cellMin, cellMax), + VofFlow::nextafter(VofFlow::nextafter(cellMax, cellMin), cellMin)); + } + + inline void cellCorrector(const VofFlow::DomainInfo& domainInfo, VofFlow::UpdateParticle& updateParticle, + vtkDataArray* vof, double eps) { + ZoneScoped; + + constexpr int searchRange = 4; + const auto startCell = domainInfo.idxToGridCoord(updateParticle.gridIdx); + + VofFlow::vec3 cellCenter_i = domainInfo.cellCenterVec3(startCell); + float minDist = std::numeric_limits::max(); + VofFlow::gridCoords_t bestCoords{-1, -1, -1}; + double bestF = -1.0; + + for (const auto& neighborCoords : VofFlow::GridRange(domainInfo.cellDims(), startCell, searchRange)) { + + VofFlow::vec3 cellCenter_n = domainInfo.cellCenterVec3(neighborCoords); + float dist = VofFlow::lengthSq(cellCenter_n - cellCenter_i); + + vtkIdType neighborIdx = domainInfo.gridCoordToIdx(neighborCoords); + double f = vof->GetComponent(neighborIdx, 0); + + if (f > eps && dist <= minDist) { + minDist = dist; + bestCoords = neighborCoords; + bestF = f; + } + } + + if (bestCoords[0] > -1) { + const auto cellMin = domainInfo.coordsVec3(bestCoords); + const auto cellMax = domainInfo.coordsVec3({bestCoords[0] + 1, bestCoords[1] + 1, bestCoords[2] + 1}); + + updateParticle.pos = clampCell(updateParticle.pos, cellMin, cellMax); + updateParticle.gridIdx = domainInfo.gridCoordToIdx(bestCoords); + updateParticle.f = bestF; + } + } + + inline void intersectPlicPlane(VofFlow::vec3& pos, const VofFlow::K::Plane_3& plane, + const VofFlow::K::Point_3& baseCorner) { + ZoneScoped; + + VofFlow::K::Point_3 pos_cgal = VofFlow::toPoint_3(pos); + if (plane.has_on_positive_side(pos_cgal)) { + VofFlow::K::Segment_3 line(baseCorner, pos_cgal); + const auto result = CGAL::intersection(plane, line); + if (result) { + VofFlow::K::Point_3 newPos; + if (const VofFlow::K::Point_3* p = boost::get(&*result)) { + newPos = *p; + } else { + const VofFlow::K::Segment_3* s = boost::get(&*result); + // This case should not be possible. Base point is always outside the plain. + // Only exception should be f = 0 cells. So just use source point of segment. + newPos = s->source(); + } + + pos = VofFlow::toVec3(newPos); + } + } + } + + inline void plicCorrector(VofFlow::UpdateParticle& updateParticle, VofFlow::CachedPlic& plicCache) { + ZoneScoped; + + const auto& plicResult = plicCache.get(updateParticle.gridIdx); + + switch (plicResult.cellClass) { + case VofFlow::CellClass::EMPTY: + case VofFlow::CellClass::FULL_PHASE1: + case VofFlow::CellClass::FULL_PHASE2: + // Should never happen. + throw std::runtime_error("PLIC Corrector called for non interface cell!"); + case VofFlow::CellClass::INTERFACE_PH1_PH2: + if (updateParticle.phaseIdx == 1) { + intersectPlicPlane(updateParticle.pos, plicResult.plic1->cell.getPlane(0), + plicResult.plic1->baseCorner); + } else { + intersectPlicPlane(updateParticle.pos, plicResult.plic1->cell.getPlane(0).opposite(), + plicResult.plic1->oppositeCorner); + } + break; + case VofFlow::CellClass::INTERFACE_PHASE1: + if (updateParticle.phaseIdx == 1) { + intersectPlicPlane(updateParticle.pos, plicResult.plic1->cell.getPlane(0), + plicResult.plic1->baseCorner); + } else { + // Should never happen. + throw std::runtime_error("PLIC Corrector called for empty phase!"); + } + break; + case VofFlow::CellClass::INTERFACE_PHASE2: + if (updateParticle.phaseIdx == 1) { + // Should never happen. + throw std::runtime_error("PLIC Corrector called for empty phase!"); + } else { + intersectPlicPlane(updateParticle.pos, plicResult.plic1->cell.getPlane(0), + plicResult.plic1->baseCorner); + } + break; + case VofFlow::CellClass::INTERFACE_ALL: + if (updateParticle.phaseIdx == 1) { + intersectPlicPlane(updateParticle.pos, plicResult.plic1->cell.getPlane(0), + plicResult.plic1->baseCorner); + } else { + // If in wrong phase, move to interface. + VofFlow::K::Point_3 pos_cgal = VofFlow::toPoint_3(updateParticle.pos); + if (plicResult.plic1->cell.getPlane(0).has_on_negative_side(pos_cgal)) { + pos_cgal = plicResult.plic1->cell.getPolygon(0).findNearestPoint(pos_cgal); + updateParticle.pos = VofFlow::toVec3(pos_cgal); + } + intersectPlicPlane(updateParticle.pos, plicResult.plic2->cell.getPlane(1), + plicResult.plic2->baseCorner); + } + break; + } + + // TODO + // Observed cases where particles are slightly out of the cell in INTERFACE_ALL cells. Looks like numerical + // instability. Clamp to cell boundary. + const auto& cell = plicResult.plic1->cell.cellCube(); + updateParticle.pos = clampCell(updateParticle.pos, VofFlow::toVec3(cell.min()), VofFlow::toVec3(cell.max())); + } +} // namespace + +void VofFlow::advectParticles(const DomainInfo& domainInfo, vtkDataArray* velocity0, vtkDataArray* velocity1, + std::vector& particles, double deltaT, int mode, int numSubSteps, const mirrorInfo_t& mirrorInfo) { + ZoneScoped; + + if (mode == 1) { +#pragma omp parallel for default(none) \ + shared(particles, domainInfo, velocity0, velocity1, deltaT, numSubSteps, mirrorInfo) + for (int64_t i = 0; i < static_cast(particles.size()); i++) { + particles[i] = rungeKutta4(particles[i], domainInfo, velocity0, velocity1, deltaT, numSubSteps, mirrorInfo); + } + } else if (mode == 2) { +#pragma omp parallel for default(none) \ + shared(particles, domainInfo, velocity0, velocity1, deltaT, numSubSteps, mirrorInfo) + for (int64_t i = 0; i < static_cast(particles.size()); i++) { + particles[i] = euler(particles[i], domainInfo, velocity0, velocity1, deltaT, numSubSteps, mirrorInfo); + } + } else { + throw std::runtime_error("Unknown Mode"); + } +} + +void VofFlow::correctParticles(const DomainInfo& domainInfo, Particles& particles, + const std::vector& oldParticlesPos, const VofData& vofData, bool neighborCorrection, bool cellCorrection, + bool plicCorrection, double eps, CachedPlic& plicCache) { + ZoneScoped; + + if (neighborCorrection == 0 && cellCorrection == 0 && plicCorrection == 0) { + return; + } + + std::unordered_map> oldCellToParticleIdMap1; + std::unordered_map> oldCellToParticleIdMap2; + if (neighborCorrection) { + for (std::size_t i = 0; i < oldParticlesPos.size(); i++) { + const auto idx = domainInfo.posToIdxEx(oldParticlesPos[i]); + if (particles.phaseIdx[i] == 1) { + oldCellToParticleIdMap1[idx].push_back(i); + } else { + oldCellToParticleIdMap2[idx].push_back(i); + } + } + } + + // Collect particles relevant for correction (particles in f < 1 cells) and store particle info. + std::vector updateParticles; + for (std::size_t p = 0; p < particles.size(); p++) { + auto pos = particles.position[p]; + bool outOfBounds = false; + auto gridCoord = domainInfo.posToGridCoord(pos); + if (!gridCoord.has_value()) { + outOfBounds = true; + // This case should be avoided, i.e. by increasing number of ghost cells. + std::cout << "Found particle out of domain!" << std::endl; // TODO + // Hotfix: pull particle back into domain. + const auto& b = domainInfo.localBounds(); + const VofFlow::vec3 cellMin{static_cast(b[0]), static_cast(b[2]), static_cast(b[4])}; + const VofFlow::vec3 cellMax{static_cast(b[1]), static_cast(b[3]), static_cast(b[5])}; + pos = clampCell(pos, cellMin, cellMax); + gridCoord = domainInfo.posToGridCoord(pos); + } + if (!gridCoord.has_value()) { + throw std::runtime_error("Invalid particle! This should not happen!"); + } + + const unsigned char phaseIdx = particles.phaseIdx[p]; + const vtkIdType idx = domainInfo.gridCoordToIdx(gridCoord.value()); + const double f = ((phaseIdx == 1) ? vofData.vof1st : vofData.vof2nd)->GetComponent(idx, 0); + if (outOfBounds || f < 1.0f - eps) { + updateParticles.push_back({p, idx, pos, phaseIdx, f}); + } + } + + // Correction + for (auto& updateParticle : updateParticles) { + const auto& phaseVof = (updateParticle.phaseIdx == 1) ? vofData.vof1st : vofData.vof2nd; + if (neighborCorrection && updateParticle.f < eps) { + const auto& phaseMap = (updateParticle.phaseIdx == 1) ? oldCellToParticleIdMap1 : oldCellToParticleIdMap2; + neighborCorrector(domainInfo, phaseMap, updateParticle, particles.position, oldParticlesPos, phaseVof, eps); + } + if (cellCorrection && updateParticle.f < eps) { + cellCorrector(domainInfo, updateParticle, phaseVof, eps); + } + if (plicCorrection && updateParticle.f > eps && updateParticle.f < 1.0 - eps) { + plicCorrector(updateParticle, plicCache); + } + } + + // Write back to particles in extra loop, because original position is needed for neighbors above. + for (auto& updateParticle : updateParticles) { + const auto idx = updateParticle.listIdx; + particles.uncertainty[idx] += length(updateParticle.pos - particles.position[idx]); + particles.position[idx] = updateParticle.pos; + } +} diff --git a/src/VofTracking/Advection.h b/src/VofTracking/Advection.h new file mode 100644 index 0000000..66819cb --- /dev/null +++ b/src/VofTracking/Advection.h @@ -0,0 +1,34 @@ +#pragma once + +#include +#include +#include + +#include + +#include "Grid/DomainInfo.h" +#include "Math/Vector.h" +#include "Misc/VofData.h" +#include "Particles.h" +#include "Plic/CachedPlic.h" + +class vtkDataArray; + +namespace VofFlow { + struct UpdateParticle { + std::size_t listIdx = 0; + vtkIdType gridIdx = 0; + vec3 pos; + unsigned char phaseIdx = 0; + double f = 0.0; + }; + + using mirrorInfo_t = std::array; + + void advectParticles(const DomainInfo& domainInfo, vtkDataArray* velocity0, vtkDataArray* velocity1, + std::vector& particles, double deltaT, int mode, int rk4NumSteps, const mirrorInfo_t& mirrorInfo); + + void correctParticles(const DomainInfo& domainInfo, Particles& particles, const std::vector& oldParticlesPos, + const VofData& vofData, bool neighborCorrection, bool cellCorrection, bool plicCorrection, double eps, + CachedPlic& plicCache); +} // namespace VofFlow diff --git a/src/VofTracking/Boundary.cpp b/src/VofTracking/Boundary.cpp new file mode 100644 index 0000000..acb83ed --- /dev/null +++ b/src/VofTracking/Boundary.cpp @@ -0,0 +1,222 @@ +#include "Boundary.h" + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Constants.h" +#include "Grid/DomainInfo.h" +#include "Grid/GridTypes.h" +#include "Misc/Profiling.h" +#include "VtkUtil/VtkUtilArray.h" + +VofFlow::BoundarySeedPoints VofFlow::exchangeBoundarySeeds(const DomainInfo& domainInfo, const SeedCoordInfo& seedInfo, + const BoundarySeedPoints& seeds, vtkMPIController* mpiController) { + ZoneScoped; + + const vtkIdType numSeeds = seeds.seedIdx->GetNumberOfTuples(); + const auto numNeighbors = domainInfo.neighbors().size(); + const int numPointsPerDim = seedInfo.numPointsPerCellDim(); + constexpr int numSharedSeeds = 2; // Number of (ghost) rows of seed points to send. + + std::vector sendSeed(numNeighbors); + + const extent_t interiorSeedExtent{ + domainInfo.localZeroExtent()[0] * numPointsPerDim + numSharedSeeds, + domainInfo.localZeroExtent()[1] * numPointsPerDim - numSharedSeeds, + domainInfo.localZeroExtent()[2] * numPointsPerDim + numSharedSeeds, + domainInfo.localZeroExtent()[3] * numPointsPerDim - numSharedSeeds, + domainInfo.localZeroExtent()[4] * numPointsPerDim + numSharedSeeds, + domainInfo.localZeroExtent()[5] * numPointsPerDim - numSharedSeeds, + }; + + std::vector neighborSeedExtents(numNeighbors); + for (std::size_t n = 0; n < numNeighbors; n++) { + const auto& neighbor = domainInfo.neighbors()[n]; + neighborSeedExtents[n] = { + neighbor.zeroExtent[0] * numPointsPerDim - numSharedSeeds, + neighbor.zeroExtent[1] * numPointsPerDim + numSharedSeeds, + neighbor.zeroExtent[2] * numPointsPerDim - numSharedSeeds, + neighbor.zeroExtent[3] * numPointsPerDim + numSharedSeeds, + neighbor.zeroExtent[4] * numPointsPerDim - numSharedSeeds, + neighbor.zeroExtent[5] * numPointsPerDim + numSharedSeeds, + }; + } + + for (vtkIdType i = 0; i < numSeeds; i++) { + // Get grid position + const auto& seedIdx = seeds.seedIdx->GetValue(i); + const auto s_coords = seedInfo.idxToSeedCoord(seedIdx); + + // Check if point is near border, to avoid loop over neighbors for interior points. + if (!Grid::isInExtent(interiorSeedExtent, s_coords)) { + for (std::size_t n = 0; n < numNeighbors; n++) { + if (Grid::isInExtent(neighborSeedExtents[n], s_coords)) { + sendSeed[n].seedIdx->InsertNextValue(seedIdx); + sendSeed[n].labels->InsertNextValue(seeds.labels->GetValue(i)); + } + } + } + } + + // Send data + const int processId = mpiController->GetLocalProcessId(); + BoundarySeedPoints result; + + for (std::size_t n = 0; n < numNeighbors; n++) { + const auto& neighbor = domainInfo.neighbors()[n]; + const auto& send = sendSeed[n]; + BoundarySeedPoints receive; + + // Pairwise send/receive between all neighbors. Process with smaller id sends first. + if (processId < neighbor.processId) { + mpiController->Send(send.seedIdx, neighbor.processId, MPITags::NEIGHBOR_SEEDS); + mpiController->Send(send.labels, neighbor.processId, MPITags::NEIGHBOR_SEEDS + 1); + mpiController->Receive(receive.seedIdx, neighbor.processId, MPITags::NEIGHBOR_SEEDS); + mpiController->Receive(receive.labels, neighbor.processId, MPITags::NEIGHBOR_SEEDS + 1); + } else { + mpiController->Receive(receive.seedIdx, neighbor.processId, MPITags::NEIGHBOR_SEEDS); + mpiController->Receive(receive.labels, neighbor.processId, MPITags::NEIGHBOR_SEEDS + 1); + mpiController->Send(send.seedIdx, neighbor.processId, MPITags::NEIGHBOR_SEEDS); + mpiController->Send(send.labels, neighbor.processId, MPITags::NEIGHBOR_SEEDS + 1); + } + + // Merge received data into single array + if (receive.seedIdx->GetNumberOfTuples() != receive.labels->GetNumberOfTuples()) { + throw std::runtime_error("Received boundary seeds have unexpected array size!"); + } + + if (receive.seedIdx->GetNumberOfTuples() > 0) { + appendArray(result.seedIdx, receive.seedIdx); + appendArray(result.labels, receive.labels); + } + } + + return result; +} + +vtkSmartPointer VofFlow::generateBoundary(const DomainInfo& domainInfo, const SeedCoordInfo& seedInfo, + const BoundarySeedPoints& seeds, const BoundarySeedPoints& neighborSeeds, int method, bool clipGhost) { + ZoneScoped; + + // Fast return if no seed points + if (seeds.seedIdx->GetNumberOfTuples() == 0) { + return vtkSmartPointer::New(); + } + + // TODO find min/max dimensions of seed points instead of using full subgrid. + + const int numPointsPerDim = seedInfo.numPointsPerCellDim(); + const auto& zeroExt = domainInfo.localZeroExtent(); + const auto& gBounds = domainInfo.globalBounds(); + + const std::array image_dim{ + (zeroExt[1] - zeroExt[0]) * numPointsPerDim + 4, + (zeroExt[3] - zeroExt[2]) * numPointsPerDim + 4, + (zeroExt[5] - zeroExt[4]) * numPointsPerDim + 4, + }; + const std::array image_min{ + zeroExt[0] * numPointsPerDim - 2, + zeroExt[2] * numPointsPerDim - 2, + zeroExt[4] * numPointsPerDim - 2, + }; + + // TODO this assumes uniform grid + const std::array spacing{ + (gBounds[1] - gBounds[0]) / static_cast(seedInfo.numPointsX()), + (gBounds[3] - gBounds[2]) / static_cast(seedInfo.numPointsY()), + (gBounds[5] - gBounds[4]) / static_cast(seedInfo.numPointsZ()), + }; + const std::array origin{ + gBounds[0] + (static_cast(image_min[0]) + 0.5) * spacing[0], + gBounds[2] + (static_cast(image_min[1]) + 0.5) * spacing[1], + gBounds[4] + (static_cast(image_min[2]) + 0.5) * spacing[2], + }; + + vtkSmartPointer data = vtkSmartPointer::New(); + data->SetDimensions(image_dim[0], image_dim[1], image_dim[2]); + data->SetOrigin(origin.data()); + data->SetSpacing(spacing.data()); + data->AllocateScalars(VTK_INT, 1); + data->GetPointData()->GetScalars()->SetName(ArrayNames::LABELS); + data->GetPointData()->GetScalars()->Fill(ErrorLabels::MAX_ERROR_LABEL); // < 0: error labels, >= 0: valid labels + + int maxLabel = -1; + + auto writeLabels = [&seedInfo, &image_min, &data, &maxLabel](const BoundarySeedPoints& s) { + for (vtkIdType i = 0; i < s.seedIdx->GetNumberOfTuples(); i++) { + const auto [s_x, s_y, s_z] = seedInfo.idxToSeedCoord(s.seedIdx->GetValue(i)); + // TODO do range check, currently relies on boundary seed exchange has no extra points + const vtkIdType idx = data->GetScalarIndex(s_x - image_min[0], s_y - image_min[1], s_z - image_min[2]); + const auto label = s.labels->GetValue(i); + data->GetPointData()->GetScalars()->SetTuple1(idx, label); + if (label > maxLabel) { + maxLabel = label; + } + } + }; + writeLabels(seeds); + writeLabels(neighborSeeds); + + vtkSmartPointer result; + if (method == 1) { + vtkNew marchingCubes; + marchingCubes->SetInputData(data); + marchingCubes->GenerateValues(maxLabel - ErrorLabels::MAX_ERROR_LABEL, ErrorLabels::MAX_ERROR_LABEL + 1, + maxLabel); + marchingCubes->Update(); + + result = marchingCubes->GetOutput(); + // Values are only shown in ParaView if the array has a name + result->GetCellData()->GetScalars()->SetName(ArrayNames::LABELS); + } else if (method == 2) { + vtkNew flyingEdges; + flyingEdges->SetInputData(data); + flyingEdges->GenerateValues(maxLabel - ErrorLabels::MAX_ERROR_LABEL, ErrorLabels::MAX_ERROR_LABEL + 1, + maxLabel); + flyingEdges->Update(); + result = flyingEdges->GetOutput(); + } else { + throw std::runtime_error("Invalid boundary method!"); + } + + // Clamp ghost cells + if (clipGhost) { + // TODO clipping a plane six times is inefficient, but using vtkBox creates holes in the surface (VTK Bug?). + vtkSmartPointer clip = vtkSmartPointer::New(); + clip->InsideOutOn(); // Exact border is kept. + vtkSmartPointer plane = vtkSmartPointer::New(); + const auto& bounds = domainInfo.localZeroBounds(); + for (int i = 0; i < 6; i++) { + const double direction = (i % 2 == 0) ? -1.0 : 1.0; + if (i / 2 == 0) { + plane->SetOrigin(bounds[i], 0.0, 0.0); + plane->SetNormal(direction, 0.0, 0.0); + } else if (i / 2 == 1) { + plane->SetOrigin(0.0, bounds[i], 0.0); + plane->SetNormal(0.0, direction, 0.0); + } else { + plane->SetOrigin(0.0, 0.0, bounds[i]); + plane->SetNormal(0.0, 0.0, direction); + } + + clip->SetInputData(result); + clip->SetClipFunction(plane); + clip->Update(); + result->ShallowCopy(clip->GetOutput()); + } + } + + return result; +} diff --git a/src/VofTracking/Boundary.h b/src/VofTracking/Boundary.h new file mode 100644 index 0000000..15f49dc --- /dev/null +++ b/src/VofTracking/Boundary.h @@ -0,0 +1,16 @@ +#pragma once + +#include +#include +#include + +#include "Grid/DomainInfo.h" +#include "Seed.h" + +namespace VofFlow { + BoundarySeedPoints exchangeBoundarySeeds(const DomainInfo& domainInfo, const SeedCoordInfo& seedInfo, + const BoundarySeedPoints& seeds, vtkMPIController* mpiController); + + vtkSmartPointer generateBoundary(const DomainInfo& domainInfo, const SeedCoordInfo& seedInfo, + const BoundarySeedPoints& seeds, const BoundarySeedPoints& neighborSeeds, int method, bool clipGhost = true); +} // namespace VofFlow diff --git a/src/VofTracking/CMakeLists.txt b/src/VofTracking/CMakeLists.txt new file mode 100644 index 0000000..5051d8c --- /dev/null +++ b/src/VofTracking/CMakeLists.txt @@ -0,0 +1,29 @@ +# Overwrite policy version inherited vom ParaView CMake. +cmake_policy(PUSH) +cmake_policy(VERSION 3.12...3.29) + +find_package(OpenMP REQUIRED) +find_package(nlohmann_json CONFIG REQUIRED) +find_package(CGAL CONFIG REQUIRED) +set(CGAL_DO_NOT_WARN_ABOUT_CMAKE_BUILD_TYPE TRUE CACHE INTERNAL "" FORCE) + +vtk_module_add_module(VofFlow::VofTracking FORCE_STATIC + CLASSES + vtkVofTracking + SOURCES + Advection.cpp + Boundary.cpp + ComponentExtraction.cpp + Particles.cpp + Seed.cpp + Uncertainty.cpp) +target_compile_features(VofTracking PUBLIC cxx_std_17) +set_target_properties(VofTracking PROPERTIES CXX_EXTENSIONS OFF) + +vtk_module_link(VofFlow::VofTracking + PRIVATE + CGAL::CGAL + nlohmann_json::nlohmann_json + OpenMP::OpenMP_CXX) + +cmake_policy(POP) diff --git a/src/VofTracking/ComponentExtraction.cpp b/src/VofTracking/ComponentExtraction.cpp new file mode 100644 index 0000000..bf4f888 --- /dev/null +++ b/src/VofTracking/ComponentExtraction.cpp @@ -0,0 +1,396 @@ +#include "ComponentExtraction.h" + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Constants.h" +#include "Grid/GridData.h" +#include "Misc/Profiling.h" +#include "VtkUtil/VtkUtilArray.h" +#include "VtkUtil/VtkUtilMPI.h" + +namespace { + struct ExtractComponentsWorker { + int nextLabelId_; + double epsilon_; + const VofFlow::DomainInfo& domainInfo_; + vtkSmartPointer labels_; + + ExtractComponentsWorker(vtkSmartPointer labels, const VofFlow::DomainInfo& domainInfo) + : nextLabelId_(0), + epsilon_(0.0), + domainInfo_(domainInfo), + labels_(std::move(labels)) {} + + inline bool isGhost(vtkIdType idx) { + return domainInfo_.ghostArray() != nullptr && domainInfo_.ghostArray()->GetValue(idx) > 0; + } + + template + void operator()(vtkAOSDataArrayTemplate* vof) { + VTK_ASSUME(vof->GetNumberOfComponents() == 1); + + vtkDataArrayAccessor> vofData(vof); + std::stack idxStack; + const vtkIdType numCells = domainInfo_.numCells(); + vtkIdType idx = 0; + + const int offsetX = 1; + const int offsetY = domainInfo_.cellDims()[0]; + const int offsetZ = domainInfo_.cellDims()[0] * domainInfo_.cellDims()[1]; + + while (idx < numCells) { + while (idx < numCells && + (vofData.Get(idx, 0) <= epsilon_ || labels_->GetValue(idx) >= 0 || isGhost(idx))) { + idx++; + } + if (idx >= numCells) { + break; + } + idxStack.push(idx); + while (!idxStack.empty()) { + const auto localIdx = idxStack.top(); + idxStack.pop(); + + if (vofData.Get(localIdx, 0) > epsilon_ && labels_->GetValue(localIdx) < 0 && !isGhost(localIdx)) { + labels_->SetValue(localIdx, nextLabelId_); + + const auto [i, j, k] = domainInfo_.idxToGridCoord(localIdx); + + if (i > 0) { + idxStack.push(localIdx - offsetX); + } + if (i < domainInfo_.cellDims()[0] - 1) { + idxStack.push(localIdx + offsetX); + } + if (j > 0) { + idxStack.push(localIdx - offsetY); + } + if (j < domainInfo_.cellDims()[1] - 1) { + idxStack.push(localIdx + offsetY); + } + if (k > 0) { + idxStack.push(localIdx - offsetZ); + } + if (k < domainInfo_.cellDims()[2] - 1) { + idxStack.push(localIdx + offsetZ); + } + } + } + nextLabelId_++; + } + } + }; +} // namespace + +VofFlow::ComponentResult VofFlow::ComponentExtractor::extractComponentsLocal(const vtkSmartPointer& vof) { + ZoneScoped; + + auto labels = + createVtkArray(ArrayNames::LABELS, 1, vof->GetNumberOfTuples(), ErrorLabels::EMPTY_COMPONENT); + + ExtractComponentsWorker worker(labels, domainInfo_); + + typedef vtkArrayDispatch::DispatchByValueType Dispatcher; + if (!Dispatcher::Execute(vof, worker)) { + throw std::runtime_error("Cannot dispatch array worker!"); + } + + return {worker.labels_, worker.nextLabelId_}; +} + +std::tuple, + std::vector> +VofFlow::ComponentExtractor::calcBorderExtents() { + ZoneScoped; + + const auto& zeroExtent = domainInfo_.localZeroExtent(); + + // Border cell rows within current zero extent. + // sendExtents are at positive dimension border, receiveExtents are at negative dimension border. + std::array sendExtents{zeroExtent, zeroExtent, zeroExtent}; + std::array receiveExtents{zeroExtent, zeroExtent, zeroExtent}; + sendExtents[0][0] = zeroExtent[1] - 1; + sendExtents[1][2] = zeroExtent[3] - 1; + sendExtents[2][4] = zeroExtent[5] - 1; + receiveExtents[0][1] = zeroExtent[0] + 1; + receiveExtents[1][3] = zeroExtent[2] + 1; + receiveExtents[2][5] = zeroExtent[4] + 1; + + // Intersect all neighbor extents to find overlapping regions + std::vector sendNeighbors; + std::vector receiveNeighbors; + + for (int i = 0; i < 3; i++) { + // Move extent by 1 to shift to neighbor zero extent + auto sendNeighborExt = sendExtents[i]; + auto receiveNeighborExt = receiveExtents[i]; + sendNeighborExt[2 * i]++; + sendNeighborExt[2 * i + 1]++; + receiveNeighborExt[2 * i]--; + receiveNeighborExt[2 * i + 1]--; + + for (const auto& neighbor : domainInfo_.neighbors()) { + auto intersection = Grid::intersection(sendNeighborExt, neighbor.zeroExtent); + if (Grid::isValidExtent(intersection)) { + // Move intersection back to current zero extent + intersection[2 * i]--; + intersection[2 * i + 1]--; + sendNeighbors.push_back({neighbor.processId, intersection}); + } + intersection = Grid::intersection(receiveNeighborExt, neighbor.zeroExtent); + if (Grid::isValidExtent(intersection)) { + // Move intersection back to current zero extent + intersection[2 * i]++; + intersection[2 * i + 1]++; + receiveNeighbors.push_back({neighbor.processId, intersection}); + } + } + } + + // Sort by process id + std::sort(sendNeighbors.begin(), sendNeighbors.end()); + std::sort(receiveNeighbors.begin(), receiveNeighbors.end()); + + return {std::move(sendNeighbors), std::move(receiveNeighbors)}; +} + +VofFlow::ComponentResult VofFlow::ComponentExtractor::extractComponents(const vtkSmartPointer& vof) { + ZoneScoped; + + auto [labels, numLabels] = extractComponentsLocal(vof); + + int numLabelsGlobal = numLabels; + + if (domainInfo_.isParallel()) { + // Idea: + // - Generate a globally unique int64 label using process id (int32) and local label (int32). + // - Calculate extent where neighbors are directly connected (top most cell layer at border, ignoring ghosts). + // - Create an ordered linear list of all cell values in these regions. + // - Send lists to neighbor, everybody sends in positive x/y/z and receives from negative x/y/z direction. + // - Compare lists and create set of connected labels. + // - Gather pairwise connections, create new globally unique labels and mapping from the combined ids + // - Send mapping to everybody, and everybody updates the label names + + const int processId = mpiController_->GetLocalProcessId(); + const int numProcesses = mpiController_->GetNumberOfProcesses(); + + const auto [sendNeighbors, receiveNeighbors] = calcBorderExtents(); + + // Collect data to send + std::vector> sendBuffers(sendNeighbors.size()); + std::vector sendBufferSizes(sendNeighbors.size()); + + for (std::size_t n = 0; n < sendNeighbors.size(); n++) { + vtkSmartPointer sendLabels = extractExtent(domainInfo_, sendNeighbors[n].extent, labels); + std::vector sendData(sendLabels->GetNumberOfTuples()); + for (vtkIdType i = 0; i < sendLabels->GetNumberOfTuples(); i++) { + sendData[i] = makeCombinedLabel(processId, sendLabels->GetValue(i)); + } + sendBufferSizes[n] = static_cast(sendData.size()); + sendBuffers[n] = std::move(sendData); + } + + // Send sizes + std::vector sendSizeRequests(sendNeighbors.size()); + for (std::size_t n = 0; n < sendNeighbors.size(); n++) { + const int sendTo = sendNeighbors[n].neighborId; + mpiController_->NoBlockSend(&sendBufferSizes[n], 1, sendTo, MPITags::COMPONENT_EXCHANGE_SIZE, + sendSizeRequests[n]); + } + + // Receive sizes + std::vector receiveBufferSizes(receiveNeighbors.size()); + std::vector receiveSizeRequests(receiveNeighbors.size()); + for (std::size_t n = 0; n < receiveNeighbors.size(); n++) { + const int receiveFrom = receiveNeighbors[n].neighborId; + mpiController_->NoBlockReceive(&receiveBufferSizes[n], 1, receiveFrom, MPITags::COMPONENT_EXCHANGE_SIZE, + receiveSizeRequests[n]); + } + + // Wait + mpiController_->WaitAll(static_cast(sendSizeRequests.size()), sendSizeRequests.data()); + mpiController_->WaitAll(static_cast(receiveSizeRequests.size()), receiveSizeRequests.data()); + + // In ParaView 5.11.1 and on systems where long is only 32-bit (Windows), the only 64-bit type for + // NoBlockSend/NoBlockReceive in vtkMPIController is vtkIdType. Therefore, just reinterpret the uint64_t ids + // as vtkIdType for sending. Keep this assert, because VTK has a build option to limit vtkIdType to 32-bit. + // ParaView 5.12 has extended the interface, so this workaround can be removed when targeting ParaView 5.12. + static_assert(sizeof(vtkIdType) == sizeof(uint64_t), "64-bit vtkIdType is required!"); + + // Send data + std::vector sendDataRequests(sendNeighbors.size()); + for (std::size_t n = 0; n < sendNeighbors.size(); n++) { + const int sendTo = sendNeighbors[n].neighborId; + mpiController_->NoBlockSend(reinterpret_cast(sendBuffers[n].data()), sendBufferSizes[n], sendTo, + MPITags::COMPONENT_EXCHANGE_DATA, sendDataRequests[n]); + } + + // Receive data + std::vector> receiveBuffers(receiveNeighbors.size()); + std::vector receiveDataRequests(receiveNeighbors.size()); + for (std::size_t n = 0; n < receiveNeighbors.size(); n++) { + const int receiveFrom = receiveNeighbors[n].neighborId; + receiveBuffers[n].resize(receiveBufferSizes[n]); + mpiController_->NoBlockReceive(reinterpret_cast(receiveBuffers[n].data()), + receiveBufferSizes[n], receiveFrom, MPITags::COMPONENT_EXCHANGE_DATA, receiveDataRequests[n]); + } + + // Wait + mpiController_->WaitAll(static_cast(sendDataRequests.size()), sendDataRequests.data()); + mpiController_->WaitAll(static_cast(receiveDataRequests.size()), receiveDataRequests.data()); + + // Compare received data + std::set sameLabels; // TODO compare to unordered set + for (std::size_t n = 0; n < receiveNeighbors.size(); n++) { + vtkSmartPointer compareLabels = extractExtent(domainInfo_, receiveNeighbors[n].extent, labels); + for (vtkIdType i = 0; i < compareLabels->GetNumberOfTuples(); i++) { + uint64_t receivedCombinedLabel = receiveBuffers[n][i]; + int localLabel = compareLabels->GetValue(i); + if (getLabelIdx(receivedCombinedLabel) >= 0 && localLabel >= 0) { + sameLabels.emplace(receivedCombinedLabel, makeCombinedLabel(processId, localLabel)); + } + } + } + + // Send all pairs to process 0 and receive local mapping + std::vector receivedMapping; + + if (processId == 0) { + std::set allLabels; + std::unordered_map> connections; + + // Process 0 + for (int l = 0; l < numLabels; l++) { + allLabels.insert(makeCombinedLabel(0, l)); + } + for (const auto& pair : sameLabels) { + connections[pair.label1].push_back(pair.label2); + connections[pair.label2].push_back(pair.label1); + } + + for (int p = 1; p < numProcesses; p++) { + int numLabels_p = 0; + mpiController_->Receive(&numLabels_p, 1, p, MPITags::COMPONENT_EXCHANGE_NUM_LABELS); + for (int l = 0; l < numLabels_p; l++) { + allLabels.insert(makeCombinedLabel(p, l)); + } + std::vector receivePairs; + ReceiveVector(mpiController_, receivePairs, p, MPITags::COMPONENT_EXCHANGE_PAIR_LIST); + for (const auto& pair : receivePairs) { + connections[pair.label1].push_back(pair.label2); + connections[pair.label2].push_back(pair.label1); + } + } + + std::vector> labelGroups; + while (!allLabels.empty()) { + std::set group; + std::queue toCheck; + toCheck.push(*allLabels.begin()); + while (!toCheck.empty()) { + uint64_t element = toCheck.front(); + toCheck.pop(); + group.insert(element); + allLabels.erase(element); + const auto it = connections.find(element); + if (it != connections.end()) { + for (uint64_t conn : it->second) { + if (group.count(conn) == 0) { + toCheck.push(conn); + } + } + connections.erase(it); // Validation + } + } + labelGroups.push_back(std::move(group)); + } + if (!connections.empty()) { + // Validation + throw std::runtime_error("Connection map not empty!"); + } + + numLabelsGlobal = static_cast(labelGroups.size()); + + // Build mapping for each process of process local label id to global label. + std::vector> mappings(numProcesses); + for (std::size_t g = 0; g < labelGroups.size(); g++) { + for (const uint64_t entry : labelGroups[g]) { + auto [label_process, label_id] = splitCombinedLabel(entry); + mappings[label_process].push_back({label_id, static_cast(g)}); + } + } + + for (int p = 1; p < numProcesses; p++) { + SendVector(mpiController_, mappings[p], p, MPITags::COMPONENT_EXCHANGE_MAPPING); + } + + receivedMapping = std::move(mappings[0]); + + } else { + mpiController_->Send(&numLabels, 1, 0, MPITags::COMPONENT_EXCHANGE_NUM_LABELS); + std::vector sendPairs(sameLabels.begin(), sameLabels.end()); + SendVector(mpiController_, sendPairs, 0, MPITags::COMPONENT_EXCHANGE_PAIR_LIST); + ReceiveVector(mpiController_, receivedMapping, 0, MPITags::COMPONENT_EXCHANGE_MAPPING); + } + + mpiController_->Broadcast(&numLabelsGlobal, 1, 0); + + // Create map from received data + std::unordered_map map; + for (const auto& m : receivedMapping) { + map[m.from] = m.to; + } + + // Update local labels to global labels. + for (vtkIdType i = 0; i < labels->GetNumberOfValues(); i++) { + int* value = labels->GetPointer(i); + if (*value >= 0) { + if (domainInfo_.isInZeroExtent(i)) { + *value = map.at(*value); + } else { + *value = ErrorLabels::BAD_COMPONENT_MAP; + } + } + } + } + + return {std::move(labels), numLabelsGlobal}; +} + +VofFlow::ComponentResult VofFlow::ComponentExtractor::checkComponents( + const vtkSmartPointer& inComponents) { + auto labels = createVtkArray(ArrayNames::LABELS, 1); + labels->ShallowCopy(inComponents); + labels->SetName(ArrayNames::LABELS); + if (labels->GetNumberOfComponents() != 1 || labels->GetNumberOfTuples() != domainInfo_.numCells()) { + throw std::runtime_error("Invalid components array!"); + } + std::array range{}; + labels->GetValueRange(range.data()); + if (range[0] < ErrorLabels::EMPTY_COMPONENT) { + throw std::runtime_error("Invalid components array with negative values!"); + } + int numLabels = range[1] + 1; + + if (domainInfo_.isParallel()) { + int numLabelsGlobal = 0; + mpiController_->AllReduce(&numLabels, &numLabelsGlobal, 1, vtkCommunicator::MAX_OP); + numLabels = numLabelsGlobal; + } + + return {labels, numLabels}; +} diff --git a/src/VofTracking/ComponentExtraction.h b/src/VofTracking/ComponentExtraction.h new file mode 100644 index 0000000..5da3a2d --- /dev/null +++ b/src/VofTracking/ComponentExtraction.h @@ -0,0 +1,90 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include "Grid/DomainInfo.h" +#include "Grid/GridTypes.h" + +namespace VofFlow { + struct ComponentResult { + vtkSmartPointer gridLabels = nullptr; + int numLabels = 0; + }; + + class ComponentExtractor { + public: + explicit ComponentExtractor(const DomainInfo& domainInfo, vtkMPIController* mpiController = nullptr) + : domainInfo_(domainInfo), + mpiController_(mpiController) {} + + ComponentResult extractComponents(const vtkSmartPointer& vof); + + ComponentResult checkComponents(const vtkSmartPointer& inComponents); + + private: + static inline uint64_t makeCombinedLabel(int processId, int label) { + return (static_cast(processId) << 32) + static_cast(label); + } + + static inline int getLabelProcessId(uint64_t l) { + return static_cast(static_cast(l >> 32)); + } + + static inline int getLabelIdx(uint64_t l) { + return static_cast(static_cast(l & 0xFFFFFFFF)); + } + + static inline std::tuple splitCombinedLabel(uint64_t l) { + return { + getLabelProcessId(l), + getLabelIdx(l), + }; + } + + struct BorderExtent { + int neighborId; + extent_t extent; + + inline bool operator<(const BorderExtent& other) const { + return neighborId < other.neighborId; + } + }; + + struct LabelPair { + uint64_t label1; + uint64_t label2; + + LabelPair() = default; + LabelPair(uint64_t l1, uint64_t l2) : label1(l1), label2(l2) { + // Label pairs should have smaller label first to avoid duplicates. + if (label2 < label1) { + std::swap(label1, label2); + } + } + + inline bool operator<(const LabelPair& other) const { + return (label1 == other.label1) ? (label2 < other.label2) : (label1 < other.label1); + } + }; + + struct MappingPair { + int from; + int to; + }; + + ComponentResult extractComponentsLocal(const vtkSmartPointer& vof); + + std::tuple, std::vector> calcBorderExtents(); + + const DomainInfo& domainInfo_; + vtkMPIController* mpiController_; + }; +} // namespace VofFlow diff --git a/src/VofTracking/Constants.h b/src/VofTracking/Constants.h new file mode 100644 index 0000000..949b47f --- /dev/null +++ b/src/VofTracking/Constants.h @@ -0,0 +1,49 @@ +#pragma once + +namespace VofFlow { + class ArrayNames { + public: + static constexpr const char* POINTS = "Points"; // Name as used by vtkPoints internally. + static constexpr const char* PHASE_IDX = "Phase"; + static constexpr const char* SEED_IDX = "SeedIdx"; + static constexpr const char* LABELS = "Labels"; + static constexpr const char* ID = "Id"; + static constexpr const char* PROCESS_ID = "ProcessId"; + static constexpr const char* CHANGED_PHASE_STEP = "ChangedPhaseStep"; + static constexpr const char* UNCERTAINTY = "Uncertainty"; + static constexpr const char* TARGET_PHASE = "TargetPhase"; + static constexpr const char* SAME_END_PHASE_RATIO = "SameEndPhaseRatio"; + static constexpr const char* STAYED_IN_PHASE_RATIO = "StayedInPhaseRatio"; + }; + + class MPITags { + public: + static constexpr int NEIGHBOR_SEEDS = 101; + static constexpr int COMPONENT_EXCHANGE_SIZE = 110; + static constexpr int COMPONENT_EXCHANGE_DATA = 111; + static constexpr int COMPONENT_EXCHANGE_NUM_LABELS = 112; + static constexpr int COMPONENT_EXCHANGE_PAIR_LIST = 113; + static constexpr int COMPONENT_EXCHANGE_MAPPING = 114; + static constexpr int PARTICLE_EXCHANGE = 120; + static constexpr int PARTICLE_TO_SEED_ID = 150; + static constexpr int PARTICLE_TO_SEED_POS = 151; + static constexpr int PARTICLE_TO_SEED_LABEL = 152; + static constexpr int PARTICLE_TO_SEED_CHANGED_PHASE_STEP = 153; + static constexpr int PARTICLE_TO_SEED_UNCERTAINTY = 154; + static constexpr int PARTICLE_TO_SEED_TARGET_PHASE = 155; + static constexpr int PARTICLE_TO_SEED_OUT_OF_BOUNDS_ID = 156; + static constexpr int PARTICLE_TO_SEED_OUT_OF_BOUNDS_CHANGED_PHASE_STEP = 157; + }; + + class ErrorLabels { + public: + static constexpr int EMPTY_COMPONENT = -1; + static constexpr int BAD_COMPONENT_MAP = -2; + static constexpr int BAD_PARTICLE = -3; + static constexpr int PARTICLE_OUT_OF_BOUNDS = -4; + static constexpr int BAD_SEED = -5; + + // Always keep last and below all other values, is used for boundary grid init. + static constexpr int MAX_ERROR_LABEL = -6; + }; +} // namespace VofFlow diff --git a/src/VofTracking/Particles.cpp b/src/VofTracking/Particles.cpp new file mode 100644 index 0000000..4ddea23 --- /dev/null +++ b/src/VofTracking/Particles.cpp @@ -0,0 +1,403 @@ +#include "Particles.h" + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "Constants.h" +#include "Grid/GridTypes.h" +#include "Misc/Profiling.h" +#include "Plic/PlicUtil.h" +#include "VtkUtil/VtkUtilArray.h" +#include "VtkUtil/VtkUtilMPI.h" + +vtkSmartPointer VofFlow::Particles::toPolyData() const { + ZoneScoped; + + vtkSmartPointer poly = vtkSmartPointer::New(); + vtkSmartPointer points = vtkSmartPointer::New(); + + points->SetData(createVtkArray(ArrayNames::POINTS, position)); + poly->SetPoints(points); + + poly->GetPointData()->AddArray(createVtkArray(ArrayNames::PHASE_IDX, phaseIdx)); + + if (!id.empty()) { + poly->GetPointData()->AddArray(createVtkArray(ArrayNames::ID, id)); + } + + if (!processId.empty()) { + poly->GetPointData()->AddArray(createVtkArray(ArrayNames::PROCESS_ID, processId)); + } + + poly->GetPointData()->AddArray(createVtkArray(ArrayNames::CHANGED_PHASE_STEP, changedPhaseStep)); + + poly->GetPointData()->AddArray(createVtkArray(ArrayNames::UNCERTAINTY, uncertainty)); + + return poly; +} + +std::unique_ptr VofFlow::initParticles(const SeedPoints& seeds, int processId) { + ZoneScoped; + + std::unique_ptr particles = std::make_unique(); + particles->position.clear(); + particles->phaseIdx.clear(); + particles->id.clear(); + particles->processId.clear(); + particles->changedPhaseStep.clear(); + particles->uncertainty.clear(); + + const vtkIdType numPoints = seeds.points->GetNumberOfTuples(); + + particles->position.resize(numPoints); + std::memcpy(particles->position.data(), seeds.points->GetVoidPointer(0), sizeof(vec3) * numPoints); + particles->phaseIdx.resize(numPoints); + std::memcpy(particles->phaseIdx.data(), seeds.phaseIdx->GetVoidPointer(0), sizeof(unsigned char) * numPoints); + // Fill id list with values from 0 to numPoints - 1. + particles->id.resize(numPoints); + std::iota(particles->id.begin(), particles->id.end(), 0); + particles->processId = std::vector(numPoints, processId); + particles->changedPhaseStep = std::vector(numPoints, -1); + particles->uncertainty = std::vector(numPoints, 0.0f); + + return particles; +} + +void VofFlow::extractOutOfBoundsParticles(Particles& particles, std::vector& particlesPosOld, + OutOfBoundsParticles& oobParticles, const DomainInfo& domainInfo, int numTimeSteps) { + ZoneScoped; + + std::size_t read_idx = 0; + std::size_t write_idx = 0; + + while (read_idx < particles.size()) { + if (domainInfo.isInGlobalBounds(particles.position[read_idx])) { + // Move elements + if (read_idx != write_idx) { + particles.position[write_idx] = particles.position[read_idx]; + particles.phaseIdx[write_idx] = particles.phaseIdx[read_idx]; + particles.id[write_idx] = particles.id[read_idx]; + particles.processId[write_idx] = particles.processId[read_idx]; + particles.changedPhaseStep[write_idx] = particles.changedPhaseStep[read_idx]; + particles.uncertainty[write_idx] = particles.uncertainty[read_idx]; + particlesPosOld[write_idx] = particlesPosOld[read_idx]; + } + read_idx++; + write_idx++; + } else { + oobParticles.id.push_back(particles.id[read_idx]); + oobParticles.processId.push_back(particles.processId[read_idx]); + const auto& s = particles.changedPhaseStep[read_idx]; + oobParticles.changedPhaseStep.push_back(s >= 0 ? s : numTimeSteps); + read_idx++; + } + } + particles.resize(write_idx); + particles.shrink_to_fit(); + particlesPosOld.resize(write_idx); + particlesPosOld.shrink_to_fit(); +} + +void VofFlow::exchangeParticles(Particles& particles, const DomainInfo& domainInfo, vtkMPIController* mpiController) { + ZoneScoped; + + const std::size_t numParticles = particles.size(); + const auto numNeighbors = domainInfo.neighbors().size(); + + std::vector sendParticles(numNeighbors); + Particles keepParticles; + // Assume most particles stay in domain + keepParticles.reserve(numParticles); + + for (std::size_t i = 0; i < numParticles; i++) { + const auto& pos = particles.position[i]; + int bound = domainInfo.isOutOfZeroBoundsDirection(pos); + if (bound < 0) { + keepParticles.position.push_back(pos); + keepParticles.phaseIdx.push_back(particles.phaseIdx[i]); + keepParticles.id.push_back(particles.id[i]); + keepParticles.processId.push_back(particles.processId[i]); + keepParticles.changedPhaseStep.push_back(particles.changedPhaseStep[i]); + keepParticles.uncertainty.push_back(particles.uncertainty[i]); + } else { + // TODO optimize neighbor search by using bound + bool found = false; + for (std::size_t n = 0; n < numNeighbors; n++) { + const auto& neighbor = domainInfo.neighbors()[n]; + if (Grid::isInBounds(neighbor.zeroBounds, {pos.x, pos.y, pos.z})) { + sendParticles[n].position.push_back(pos); + sendParticles[n].phaseIdx.push_back(particles.phaseIdx[i]); + sendParticles[n].id.push_back(particles.id[i]); + sendParticles[n].processId.push_back(particles.processId[i]); + sendParticles[n].changedPhaseStep.push_back(particles.changedPhaseStep[i]); + sendParticles[n].uncertainty.push_back(particles.uncertainty[i]); + found = true; + break; + } + } + if (!found) { + throw std::runtime_error("Particle left domain from neighbors!"); + } + } + } + particles = std::move(keepParticles); + + // Send / receive + const int processId = mpiController->GetLocalProcessId(); + std::vector receiveParticles(numNeighbors); + + for (std::size_t n = 0; n < numNeighbors; n++) { + const auto& neighbor = domainInfo.neighbors()[n]; + const auto& send = sendParticles[n]; + auto& receive = receiveParticles[n]; + + // Pairwise send/receive between all neighbors. Process with smaller id sends first. + if (processId < neighbor.processId) { + SendVector(mpiController, send.position, neighbor.processId, MPITags::PARTICLE_EXCHANGE); + SendVector(mpiController, send.phaseIdx, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 1); + SendVector(mpiController, send.id, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 2); + SendVector(mpiController, send.processId, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 3); + SendVector(mpiController, send.changedPhaseStep, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 4); + SendVector(mpiController, send.uncertainty, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 5); + ReceiveVector(mpiController, receive.position, neighbor.processId, MPITags::PARTICLE_EXCHANGE); + ReceiveVector(mpiController, receive.phaseIdx, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 1); + ReceiveVector(mpiController, receive.id, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 2); + ReceiveVector(mpiController, receive.processId, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 3); + ReceiveVector(mpiController, receive.changedPhaseStep, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 4); + ReceiveVector(mpiController, receive.uncertainty, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 5); + } else { + ReceiveVector(mpiController, receive.position, neighbor.processId, MPITags::PARTICLE_EXCHANGE); + ReceiveVector(mpiController, receive.phaseIdx, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 1); + ReceiveVector(mpiController, receive.id, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 2); + ReceiveVector(mpiController, receive.processId, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 3); + ReceiveVector(mpiController, receive.changedPhaseStep, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 4); + ReceiveVector(mpiController, receive.uncertainty, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 5); + SendVector(mpiController, send.position, neighbor.processId, MPITags::PARTICLE_EXCHANGE); + SendVector(mpiController, send.phaseIdx, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 1); + SendVector(mpiController, send.id, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 2); + SendVector(mpiController, send.processId, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 3); + SendVector(mpiController, send.changedPhaseStep, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 4); + SendVector(mpiController, send.uncertainty, neighbor.processId, MPITags::PARTICLE_EXCHANGE + 5); + } + } + + for (std::size_t n = 0; n < numNeighbors; n++) { + const auto& receive = receiveParticles[n]; + if (!receive.hasValidSize()) { + throw std::runtime_error("Received bad particle array sizes"); + } + + particles.position.insert(particles.position.end(), receive.position.begin(), receive.position.end()); + particles.phaseIdx.insert(particles.phaseIdx.end(), receive.phaseIdx.begin(), receive.phaseIdx.end()); + particles.id.insert(particles.id.end(), receive.id.begin(), receive.id.end()); + particles.processId.insert(particles.processId.end(), receive.processId.begin(), receive.processId.end()); + particles.changedPhaseStep.insert(particles.changedPhaseStep.end(), receive.changedPhaseStep.begin(), + receive.changedPhaseStep.end()); + particles.uncertainty.insert(particles.uncertainty.end(), receive.uncertainty.begin(), + receive.uncertainty.end()); + } + particles.shrink_to_fit(); +} + +std::vector VofFlow::labelAdvectedParticles(const DomainInfo& domainInfo, const ComponentResult& components1, + const ComponentResult& components2, const Particles& particles, vtkImplicitFunction* cutFunction) { + ZoneScoped; + + std::vector labels(particles.size()); + + for (std::size_t i = 0; i < particles.size(); ++i) { + const auto& pos = particles.position[i]; + const auto phase = particles.phaseIdx[i]; + const auto& cell = domainInfo.posToGridCoord(pos); + int label = ErrorLabels::BAD_PARTICLE; + if (cell.has_value()) { + if (phase == 1) { + label = components1.gridLabels->GetValue(domainInfo.gridCoordToIdx(cell.value())); + } else { + if (components2.gridLabels == nullptr) { + throw std::runtime_error("Invalid phase!"); + } + label = components2.gridLabels->GetValue(domainInfo.gridCoordToIdx(cell.value())); + if (label >= 0) { + // Offset to make labels in second phase unique. + label += components1.numLabels; + } + } + if (label >= 0 && cutFunction != nullptr && cutFunction->FunctionValue(pos.x, pos.y, pos.z) < 0) { + // Offset to make cut regions unique labels. + label += (components1.numLabels + components2.numLabels); + } + } + labels[i] = label; + } + + return labels; +} + +void VofFlow::checkPhaseChanged(Particles& particles, CachedPlic& plicCache, int numTimeSteps) { + ZoneScoped; + + for (std::size_t p = 0; p < particles.size(); p++) { + // Ignore particles which already changed phase + if (particles.changedPhaseStep[p] >= 0) { + continue; + } + const auto& pos = particles.position[p]; + auto plicResult = plicCache.get(pos); + if (particles.phaseIdx[p] != getPhaseIdx(plicResult, pos)) { + particles.changedPhaseStep[p] = numTimeSteps; + } + } +} + +std::vector VofFlow::getParticlePhase(CachedPlic& plicCache, const Particles& particles) { + ZoneScoped; + + std::vector phaseIdxResult(particles.size()); + + for (std::size_t p = 0; p < particles.size(); p++) { + const auto& pos = particles.position[p]; + auto plicResult = plicCache.get(pos); + + phaseIdxResult[p] = getPhaseIdx(plicResult, pos); + } + + return phaseIdxResult; +} + +VofFlow::SeedResult VofFlow::transferParticleDataToSeeds(const SeedPoints& seeds, const Particles& particles, + const std::vector& labels, const std::vector& targetPhases, + const OutOfBoundsParticles& oobParticles, vtkMPIController* mpiController) { + ZoneScoped; + + SeedResult result(seeds.points->GetNumberOfTuples()); + + if (mpiController != nullptr) { + const int processId = mpiController->GetLocalProcessId(); + const int numProcesses = mpiController->GetNumberOfProcesses(); + + std::vector> sendId(numProcesses); + std::vector> sendPosition(numProcesses); + std::vector> sendLabel(numProcesses); + std::vector> sendChangedPhaseStep(numProcesses); + std::vector> sendUncertainty(numProcesses); + std::vector> sendTargetPhase(numProcesses); + std::vector> sendOutOfBoundsId(numProcesses); + std::vector> sendOutOfBoundsChangedPhaseStep(numProcesses); + + for (std::size_t i = 0; i < particles.size(); i++) { + const auto& sendProcess = particles.processId[i]; + sendId[sendProcess].push_back(particles.id[i]); + sendPosition[sendProcess].push_back(particles.position[i]); + sendLabel[sendProcess].push_back(labels[i]); + sendChangedPhaseStep[sendProcess].push_back(particles.changedPhaseStep[i]); + sendUncertainty[sendProcess].push_back(particles.uncertainty[i]); + sendTargetPhase[sendProcess].push_back(targetPhases[i]); + } + for (std::size_t i = 0; i < oobParticles.id.size(); i++) { + const auto& sendProcess = oobParticles.processId[i]; + sendOutOfBoundsId[sendProcess].push_back(oobParticles.id[i]); + sendOutOfBoundsChangedPhaseStep[sendProcess].push_back(oobParticles.changedPhaseStep[i]); + } + + std::vector> receiveId(numProcesses); + std::vector> receivePosition(numProcesses); + std::vector> receiveLabel(numProcesses); + std::vector> receiveChangedPhaseStep(numProcesses); + std::vector> receiveUncertainty(numProcesses); + std::vector> receiveTargetPhase(numProcesses); + std::vector> receiveOutOfBoundsId(numProcesses); + std::vector> receiveOutOfBoundsChangedPhaseStep(numProcesses); + + for (int i = 0; i < numProcesses - 1; i++) { + // Send all to all, everyone starts with process id +1 neighbor. + int sendTo = (processId + i + 1) % numProcesses; + int receiveFrom = (processId - i - 1 + numProcesses) % numProcesses; + + SendVector(mpiController, sendId[sendTo], sendTo, MPITags::PARTICLE_TO_SEED_ID); + SendVector(mpiController, sendPosition[sendTo], sendTo, MPITags::PARTICLE_TO_SEED_POS); + SendVector(mpiController, sendLabel[sendTo], sendTo, MPITags::PARTICLE_TO_SEED_LABEL); + SendVector(mpiController, sendChangedPhaseStep[sendTo], sendTo, + MPITags::PARTICLE_TO_SEED_CHANGED_PHASE_STEP); + SendVector(mpiController, sendUncertainty[sendTo], sendTo, MPITags::PARTICLE_TO_SEED_UNCERTAINTY); + SendVector(mpiController, sendTargetPhase[sendTo], sendTo, MPITags::PARTICLE_TO_SEED_TARGET_PHASE); + SendVector(mpiController, sendOutOfBoundsId[sendTo], sendTo, MPITags::PARTICLE_TO_SEED_OUT_OF_BOUNDS_ID); + SendVector(mpiController, sendOutOfBoundsChangedPhaseStep[sendTo], sendTo, + MPITags::PARTICLE_TO_SEED_OUT_OF_BOUNDS_CHANGED_PHASE_STEP); + ReceiveVector(mpiController, receiveId[receiveFrom], receiveFrom, MPITags::PARTICLE_TO_SEED_ID); + ReceiveVector(mpiController, receivePosition[receiveFrom], receiveFrom, MPITags::PARTICLE_TO_SEED_POS); + ReceiveVector(mpiController, receiveLabel[receiveFrom], receiveFrom, MPITags::PARTICLE_TO_SEED_LABEL); + ReceiveVector(mpiController, receiveChangedPhaseStep[receiveFrom], receiveFrom, + MPITags::PARTICLE_TO_SEED_CHANGED_PHASE_STEP); + ReceiveVector(mpiController, receiveUncertainty[receiveFrom], receiveFrom, + MPITags::PARTICLE_TO_SEED_UNCERTAINTY); + ReceiveVector(mpiController, receiveTargetPhase[receiveFrom], receiveFrom, + MPITags::PARTICLE_TO_SEED_TARGET_PHASE); + ReceiveVector(mpiController, receiveOutOfBoundsId[receiveFrom], receiveFrom, + MPITags::PARTICLE_TO_SEED_OUT_OF_BOUNDS_ID); + ReceiveVector(mpiController, receiveOutOfBoundsChangedPhaseStep[receiveFrom], receiveFrom, + MPITags::PARTICLE_TO_SEED_OUT_OF_BOUNDS_CHANGED_PHASE_STEP); + } + + // Copy self + std::swap(receiveId[processId], sendId[processId]); + std::swap(receivePosition[processId], sendPosition[processId]); + std::swap(receiveLabel[processId], sendLabel[processId]); + std::swap(receiveChangedPhaseStep[processId], sendChangedPhaseStep[processId]); + std::swap(receiveUncertainty[processId], sendUncertainty[processId]); + std::swap(receiveTargetPhase[processId], sendTargetPhase[processId]); + std::swap(receiveOutOfBoundsId[processId], sendOutOfBoundsId[processId]); + std::swap(receiveOutOfBoundsChangedPhaseStep[processId], sendOutOfBoundsChangedPhaseStep[processId]); + + for (int p = 0; p < numProcesses; p++) { + for (std::size_t i = 0; i < receiveId[p].size(); i++) { + int idx = receiveId[p][i]; + result.positions->SetTypedTuple(idx, receivePosition[p][i].data()); + result.labels->SetTypedComponent(idx, 0, receiveLabel[p][i]); + result.changedPhaseStep->SetTypedComponent(idx, 0, receiveChangedPhaseStep[p][i]); + result.uncertainty->SetTypedComponent(idx, 0, receiveUncertainty[p][i]); + result.targetPhase->SetTypedComponent(idx, 0, receiveTargetPhase[p][i]); + } + for (std::size_t i = 0; i < receiveOutOfBoundsId[p].size(); i++) { + constexpr std::array zero{0.0f, 0.0f, 0.0f}; + int idx = receiveOutOfBoundsId[p][i]; + result.positions->SetTypedTuple(idx, zero.data()); + result.labels->SetTypedComponent(idx, 0, ErrorLabels::PARTICLE_OUT_OF_BOUNDS); + result.changedPhaseStep->SetTypedComponent(idx, 0, receiveOutOfBoundsChangedPhaseStep[p][i]); + result.uncertainty->SetTypedComponent(idx, 0, -1.0); + result.targetPhase->SetTypedComponent(idx, 0, 0); + } + } + } else { + // As particles may go out of bounds, need to check one by one. + for (std::size_t p = 0; p < particles.size(); p++) { + const auto& idx = particles.id[p]; + result.positions->SetTypedTuple(idx, particles.position[p].data()); + result.labels->SetTypedComponent(idx, 0, labels[p]); + result.changedPhaseStep->SetTypedComponent(idx, 0, particles.changedPhaseStep[p]); + result.uncertainty->SetTypedComponent(idx, 0, particles.uncertainty[p]); + result.targetPhase->SetTypedComponent(idx, 0, targetPhases[p]); + } + for (std::size_t i = 0; i < oobParticles.id.size(); i++) { + const auto idx = oobParticles.id[i]; + constexpr std::array zero{0.0f, 0.0f, 0.0f}; + result.positions->SetTypedTuple(idx, zero.data()); + result.labels->SetTypedComponent(idx, 0, ErrorLabels::PARTICLE_OUT_OF_BOUNDS); + result.changedPhaseStep->SetTypedComponent(idx, 0, oobParticles.changedPhaseStep[i]); + result.uncertainty->SetTypedComponent(idx, 0, -1.0); + result.targetPhase->SetTypedComponent(idx, 0, 0); + } + } + + return result; +} diff --git a/src/VofTracking/Particles.h b/src/VofTracking/Particles.h new file mode 100644 index 0000000..0499cc2 --- /dev/null +++ b/src/VofTracking/Particles.h @@ -0,0 +1,107 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include +#include + +#include "ComponentExtraction.h" +#include "Grid/DomainInfo.h" +#include "Math/Vector.h" +#include "Plic/CachedPlic.h" +#include "Seed.h" + +namespace VofFlow { + struct Particles { + // With some quick empiric testing the MPI communication overhead of VTK arrays is higher than the conversion + // from std::vector to VTK arrays for a larger number of MPI processes (order of magnitude ~64). Using VTK + // arrays directly is faster for single processes or MPI processes up to ~16, but we want to optimize for a + // larger number of processes here. Therefore, use std::vector structures for particle data. + + std::vector position; + std::vector phaseIdx; + std::vector id; + std::vector processId; + std::vector changedPhaseStep; + std::vector uncertainty; + + Particles() = default; + + [[nodiscard]] vtkSmartPointer toPolyData() const; + + bool hasValidSize() const { + const auto s = position.size(); + return s == phaseIdx.size() && s == id.size() && s == processId.size() && s == changedPhaseStep.size() && + s == uncertainty.size(); + } + + [[nodiscard]] std::size_t size() const { + // assume all arrays are same + return position.size(); + } + + void clear() { + position.clear(); + phaseIdx.clear(); + id.clear(); + processId.clear(); + changedPhaseStep.clear(); + uncertainty.clear(); + } + + void reserve(std::size_t size) { + position.reserve(size); + phaseIdx.reserve(size); + id.reserve(size); + processId.reserve(size); + changedPhaseStep.reserve(size); + uncertainty.reserve(size); + } + + void resize(std::size_t size) { + position.resize(size); + phaseIdx.resize(size); + id.resize(size); + processId.resize(size); + changedPhaseStep.resize(size); + uncertainty.resize(size); + } + + void shrink_to_fit() { + position.shrink_to_fit(); + phaseIdx.shrink_to_fit(); + id.shrink_to_fit(); + processId.shrink_to_fit(); + changedPhaseStep.shrink_to_fit(); + uncertainty.shrink_to_fit(); + } + }; + + struct OutOfBoundsParticles { + std::vector id; + std::vector processId; + std::vector changedPhaseStep; + }; + + std::unique_ptr initParticles(const SeedPoints& seeds, int processId); + + void extractOutOfBoundsParticles(Particles& particles, std::vector& particlesPosOld, + OutOfBoundsParticles& oobParticles, const DomainInfo& domainInfo, int numTimeSteps); + + void exchangeParticles(Particles& particles, const DomainInfo& domainInfo, vtkMPIController* mpiController); + + std::vector labelAdvectedParticles(const DomainInfo& domainInfo, const ComponentResult& components1, + const ComponentResult& components2, const Particles& particles, vtkImplicitFunction* cutFunction); + + void checkPhaseChanged(Particles& particles, CachedPlic& plicCache, int numTimeSteps); + + std::vector getParticlePhase(CachedPlic& plicCache, const Particles& particles); + + SeedResult transferParticleDataToSeeds(const SeedPoints& seeds, const Particles& particles, + const std::vector& labels, const std::vector& targetPhases, + const OutOfBoundsParticles& oobParticles, vtkMPIController* mpiController); +} // namespace VofFlow diff --git a/src/VofTracking/Seed.cpp b/src/VofTracking/Seed.cpp new file mode 100644 index 0000000..c4b5996 --- /dev/null +++ b/src/VofTracking/Seed.cpp @@ -0,0 +1,133 @@ +#include "Seed.h" + +#include +#include +#include + +#include +#include +#include + +#include "Constants.h" +#include "Grid/GridIterator.h" +#include "Math/Vector.h" +#include "Misc/Profiling.h" +#include "Plic/PlicUtil.h" +#include "VtkUtil/VtkUtilArray.h" + +VofFlow::SeedCoordInfo::SeedCoordInfo(const DomainInfo& domainInfo, int refinement) + : domainInfo_(domainInfo), + refinement_(refinement) { + numPointsPerCellDim_ = calcNumPointsPerCellDim(refinement_); + const auto globalDims = Grid::extentDimensions(domainInfo_.globalExtent()); + numPointsX_ = static_cast(globalDims[0]) * static_cast(numPointsPerCellDim_); + numPointsY_ = static_cast(globalDims[1]) * static_cast(numPointsPerCellDim_); + numPointsZ_ = static_cast(globalDims[2]) * static_cast(numPointsPerCellDim_); +} + +VofFlow::SeedPoints::SeedPoints() { + points = createVtkArray(ArrayNames::POINTS, 3); + phaseIdx = createVtkArray(ArrayNames::PHASE_IDX, 1); + seedIdx = createVtkArray(ArrayNames::SEED_IDX, 1); +} + +VofFlow::SeedResult::SeedResult(vtkIdType numTuples) { + positions = createVtkArray(ArrayNames::POINTS, 3, numTuples, 0.0); + labels = createVtkArray(ArrayNames::LABELS, 1, numTuples, ErrorLabels::BAD_SEED); + changedPhaseStep = createVtkArray(ArrayNames::CHANGED_PHASE_STEP, 1, numTuples, -1); + uncertainty = createVtkArray(ArrayNames::UNCERTAINTY, 1, numTuples, -1); + targetPhase = createVtkArray(ArrayNames::TARGET_PHASE, 1, numTuples, 0); +} + +VofFlow::BoundarySeedPoints::BoundarySeedPoints() { + seedIdx = createVtkArray(ArrayNames::SEED_IDX, 1); + labels = createVtkArray(ArrayNames::LABELS, 1); +} + +VofFlow::BoundarySeedPoints::BoundarySeedPoints(vtkSmartPointer i, vtkSmartPointer l) + : seedIdx(std::move(i)), + labels(std::move(l)) { + if (seedIdx == nullptr || labels == nullptr || seedIdx->GetNumberOfTuples() != labels->GetNumberOfTuples()) { + throw std::runtime_error("Invalid BoundarySeedPoints arrays!"); + } +} + +vtkSmartPointer VofFlow::makePolyData(const SeedPoints& seedPoints, const SeedResult& seedResult) { + ZoneScoped; + + // Assume that SeedPoints and SeedResult have consistent array sizes within. + if (seedPoints.points->GetNumberOfTuples() != seedResult.positions->GetNumberOfTuples()) { + throw std::runtime_error("Invalid array sizes for poly data!"); + } + + vtkSmartPointer data = vtkSmartPointer::New(); + vtkSmartPointer data_points = vtkSmartPointer::New(); + data_points->SetData(seedPoints.points); + data->SetPoints(data_points); + data->GetPointData()->AddArray(seedPoints.phaseIdx); + data->GetPointData()->AddArray(seedPoints.seedIdx); + // ignore seedResult.positions + data->GetPointData()->AddArray(seedResult.labels); + data->GetPointData()->AddArray(seedResult.changedPhaseStep); + data->GetPointData()->AddArray(seedResult.uncertainty); + data->GetPointData()->AddArray(seedResult.targetPhase); + return data; +} + +std::unique_ptr VofFlow::generateSeedPoints(const DomainInfo& domainInfo, + const SeedCoordInfo& seedInfo, const VofData& vofData, double eps, int numIterations) { + ZoneScoped; + + // Make list of relative point positions within each cell along one dimension in [0, 1] range. + const int numPointsPerDim = seedInfo.numPointsPerCellDim(); + std::vector offset(numPointsPerDim); + for (int i = 0; i < numPointsPerDim; i++) { + offset[i] = (static_cast(i) + 0.5f) / static_cast(numPointsPerDim); + } + + auto seedPoints = std::make_unique(); + + // Loop over zero extent cells + const auto& extLocal = domainInfo.localExtent(); + const auto& extZero = domainInfo.localZeroExtent(); + for (const gridCoords_t& e_coords : GridRange(extZero)) { + // Grid coords + const int g_x = e_coords[0] - extLocal[0]; + const int g_y = e_coords[1] - extLocal[2]; + const int g_z = e_coords[2] - extLocal[4]; + + const auto gridIdx = domainInfo.gridCoordToIdx(g_x, g_y, g_z); + + const auto& plicResult = calcPlicOrPlic3CellClass(domainInfo, {g_x, g_y, g_z}, vofData, eps, numIterations); + + if (plicResult.cellClass == CellClass::EMPTY) { + continue; + } + + const auto cellStart = domainInfo.coordsVec3(g_x, g_y, g_z); + const auto cellSize = domainInfo.cellSizeVec3(g_x, g_y, g_z); + for (int x = 0; x < numPointsPerDim; x++) { + for (int y = 0; y < numPointsPerDim; y++) { + for (int z = 0; z < numPointsPerDim; z++) { + const vec3 pos = cellStart + cellSize * vec3(offset[x], offset[y], offset[z]); + const unsigned char phase = getPhaseIdx(plicResult, pos); + if (phase == 0) { + continue; + } + + // Seed coords + const int s_x = e_coords[0] * numPointsPerDim + x; + const int s_y = e_coords[1] * numPointsPerDim + y; + const int s_z = e_coords[2] * numPointsPerDim + z; + + seedPoints->gridIdxToSeedsMap[gridIdx].push_back(seedPoints->points->GetNumberOfTuples()); + seedPoints->points->InsertNextTypedTuple(pos.data()); + seedPoints->phaseIdx->InsertNextValue(phase); + seedPoints->seedIdx->InsertNextValue(seedInfo.seedCoordToIdx(s_x, s_y, s_z)); + } + } + } + } + + return seedPoints; +} diff --git a/src/VofTracking/Seed.h b/src/VofTracking/Seed.h new file mode 100644 index 0000000..fa526f7 --- /dev/null +++ b/src/VofTracking/Seed.h @@ -0,0 +1,99 @@ +#pragma once + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include "Grid/DomainInfo.h" +#include "Grid/GridTypes.h" +#include "Misc/VofData.h" + +namespace VofFlow { + class SeedCoordInfo { + public: + SeedCoordInfo(const DomainInfo& domainInfo, int refinement); + + [[nodiscard]] int numPointsPerCellDim() const { + return numPointsPerCellDim_; + } + + [[nodiscard]] vtkIdType numPointsX() const { + return numPointsX_; + } + + [[nodiscard]] vtkIdType numPointsY() const { + return numPointsY_; + } + + [[nodiscard]] vtkIdType numPointsZ() const { + return numPointsZ_; + } + + [[nodiscard]] inline vtkIdType seedCoordToIdx(int s_x, int s_y, int s_z) const { + return static_cast(s_x) + numPointsX_ * static_cast(s_y) + + numPointsX_ * numPointsY_ * static_cast(s_z); + } + + [[nodiscard]] inline gridCoords_t idxToSeedCoord(vtkIdType idx) const { + const int x = static_cast(idx % numPointsX_); + const int y = static_cast((idx / numPointsX_) % numPointsY_); + const int z = static_cast(idx / (numPointsX_ * numPointsY_)); + return {x, y, z}; + } + + private: + static inline int calcNumPointsPerCellDim(int refinement) { + return std::max(0, refinement) + 1; + } + + const DomainInfo& domainInfo_; + int refinement_; + int numPointsPerCellDim_; + vtkIdType numPointsX_; + vtkIdType numPointsY_; + vtkIdType numPointsZ_; + }; + + struct SeedPoints { + vtkSmartPointer points; + vtkSmartPointer phaseIdx; + vtkSmartPointer seedIdx; + + SeedPoints(); + + std::unordered_map> gridIdxToSeedsMap; + }; + + struct SeedResult { + vtkSmartPointer positions; + vtkSmartPointer labels; + vtkSmartPointer changedPhaseStep; + vtkSmartPointer uncertainty; + vtkSmartPointer targetPhase; + + explicit SeedResult(vtkIdType numTuples = 0); + }; + + struct BoundarySeedPoints { + vtkSmartPointer seedIdx; + vtkSmartPointer labels; + + BoundarySeedPoints(); + + BoundarySeedPoints(vtkSmartPointer i, vtkSmartPointer l); + }; + + vtkSmartPointer makePolyData(const SeedPoints& seedPoints, const SeedResult& seedResult); + + std::unique_ptr generateSeedPoints(const DomainInfo& domainInfo, const SeedCoordInfo& seedInfo, + const VofData& vofData, double eps, int numIterations); +} // namespace VofFlow diff --git a/src/VofTracking/TimeMeasure.h b/src/VofTracking/TimeMeasure.h new file mode 100644 index 0000000..6947dd1 --- /dev/null +++ b/src/VofTracking/TimeMeasure.h @@ -0,0 +1,167 @@ +#pragma once + +#ifndef MEASURE_TIME + +#define TIME_MEASURE_RESET(C) +#define TIME_MEASURE_START(N) +#define TIME_MEASURE_END(N) +#define TIME_MEASURE_SYNC_START(N) +#define TIME_MEASURE_SYNC_END(N) +#define TIME_MEASURE_BARRIER() +#define TIME_MEASURE_JSON(J) + +#else + +#define TIME_MEASURE_RESET(C) TimeMeasure::getInstance().reset(C) +#define TIME_MEASURE_START(N) TimeMeasure::getInstance().start(N) +#define TIME_MEASURE_END(N) TimeMeasure::getInstance().end(N) +#define TIME_MEASURE_SYNC_START(N) TimeMeasure::getInstance().syncStart(N) +#define TIME_MEASURE_SYNC_END(N) TimeMeasure::getInstance().syncEnd(N) +#define TIME_MEASURE_BARRIER() TimeMeasure::getInstance().barrier() +#define TIME_MEASURE_JSON(J) TimeMeasure::getInstance().toJson(J) + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +class TimeMeasure { +public: + using json = nlohmann::json; + + enum class Name { + All = 0, + Seeding, + Advection, + Advection_Init, + Advection_Advection, + Advection_OOB, + Advection_Corrector, + Advection_PhaseChange, + Advection_Exchange, + Components, + ParticleLabeling, + SeedLabeling, + Boundary, + Output, + }; + + static inline TimeMeasure& getInstance() { + if (instance_ == nullptr) { + instance_ = std::unique_ptr(new TimeMeasure()); + } + return *instance_; + } + + inline void reset(vtkMPIController* mpiController = nullptr) { + mpiController_ = mpiController; + // Store a high precision time stamp for measuring nanoseconds relative to this and a system clock to print + // actual date and time of measurement start. The actual date and time is meant as metadata of the measurement + // and should not be combined with the relative nanoseconds. Therefore, offset between them is ignored here. + barrier(); + resetTime_ = clock_type::now(); + resetTimeSystem_ = std::chrono::system_clock::now(); + start_.clear(); + end_.clear(); + } + + inline void barrier() { + if (mpiController_ != nullptr && mpiController_->GetCommunicator() != nullptr) { + mpiController_->Barrier(); + } + } + + inline void start(Name name) { + start_[name].push_back(clock_type::now()); + } + + inline void end(Name name) { + end_[name].push_back(clock_type::now()); + } + + inline void syncStart(Name name) { + barrier(); + start(name); + } + + inline void syncEnd(Name name) { + barrier(); + end(name); + } + + void toJson(json& j) { + std::stringstream s; + auto time = std::chrono::system_clock::to_time_t(resetTimeSystem_); + s << std::put_time(std::localtime(&time), "%FT%T%z"); + j["ResetTime"] = s.str(); + + json timings; + timings["All"] = getTimingsJson(Name::All); + timings["Seeding"] = getTimingsJson(Name::Seeding); + timings["Advection"] = getTimingsJson(Name::Advection); + timings["AdvectionInit"] = getTimingsJson(Name::Advection_Init); + timings["AdvectionAdvection"] = getTimingsJson(Name::Advection_Advection); + timings["AdvectionOOB"] = getTimingsJson(Name::Advection_OOB); + timings["AdvectionCorrector"] = getTimingsJson(Name::Advection_Corrector); + timings["AdvectionPhaseChange"] = getTimingsJson(Name::Advection_PhaseChange); + timings["AdvectionExchange"] = getTimingsJson(Name::Advection_Exchange); + timings["Components"] = getTimingsJson(Name::Components); + timings["ParticleLabeling"] = getTimingsJson(Name::ParticleLabeling); + timings["SeedLabeling"] = getTimingsJson(Name::SeedLabeling); + timings["Boundary"] = getTimingsJson(Name::Boundary); + timings["Output"] = getTimingsJson(Name::Output); + j["Timings"] = timings; + } + +private: + typedef std::chrono::high_resolution_clock clock_type; + typedef std::unordered_map> map_type; + + TimeMeasure() : mpiController_(nullptr) { + reset(); + } + + inline std::vector nanoSecondsSinceReset(const std::vector& time) const { + std::vector result(time.size()); + for (std::size_t i = 0; i < time.size(); i++) { + result[i] = std::chrono::duration_cast(time[i] - resetTime_).count(); + } + return result; + } + + inline json getTimingsJson(Name name) const { + auto start_it = start_.find(name); + auto end_it = end_.find(name); + if (start_it == start_.end() || end_it == end_.end()) { + return {}; + } + auto start = start_it->second; + auto end = end_it->second; + if (start.size() != end.size()) { + throw std::runtime_error("Invalid time measure!"); + } + return { + {"start", nanoSecondsSinceReset(start)}, + {"end", nanoSecondsSinceReset(end)}, + }; + } + + static inline std::unique_ptr instance_ = nullptr; + vtkMPIController* mpiController_; + clock_type::time_point resetTime_; + std::chrono::system_clock::time_point resetTimeSystem_; + map_type start_; + map_type end_; +}; + +#endif diff --git a/src/VofTracking/Uncertainty.cpp b/src/VofTracking/Uncertainty.cpp new file mode 100644 index 0000000..6a79bbe --- /dev/null +++ b/src/VofTracking/Uncertainty.cpp @@ -0,0 +1,46 @@ +#include "Uncertainty.h" + +#include +#include + +#include + +#include "Constants.h" +#include "Grid/GridIterator.h" +#include "Misc/Profiling.h" +#include "VtkUtil/VtkUtilArray.h" + +VofFlow::UncertaintyArrays VofFlow::makeUncertaintyArrays(const DomainInfo& domainInfo, const SeedPoints& seeds, + const SeedResult& seedResult) { + ZoneScoped; + + UncertaintyArrays result; + result.sameEndPhaseRatio = + createVtkArray(ArrayNames::SAME_END_PHASE_RATIO, 1, domainInfo.numCells(), -1.0); + result.stayedInPhaseRatio = + createVtkArray(ArrayNames::STAYED_IN_PHASE_RATIO, 1, domainInfo.numCells(), -1.0); + + for (const auto g_coords : GridRange({0, 0, 0}, domainInfo.cellDims())) { + const auto gridIdx = domainInfo.gridCoordToIdx(g_coords); + auto it = seeds.gridIdxToSeedsMap.find(gridIdx); + if (it == seeds.gridIdxToSeedsMap.end()) { + continue; + } + int numTotal = 0; + int numSameEndPhase = 0; + int numStayedInPhase = 0; + for (const auto& seedListIdx : it->second) { + numTotal++; + if (seeds.phaseIdx->GetValue(seedListIdx) == seedResult.targetPhase->GetValue(seedListIdx)) { + numSameEndPhase++; + } + if (seedResult.changedPhaseStep->GetValue(seedListIdx) < 0) { + numStayedInPhase++; + } + } + result.sameEndPhaseRatio->SetValue(gridIdx, static_cast(numSameEndPhase) / static_cast(numTotal)); + result.stayedInPhaseRatio->SetValue(gridIdx, + static_cast(numStayedInPhase) / static_cast(numTotal)); + } + return result; +} diff --git a/src/VofTracking/Uncertainty.h b/src/VofTracking/Uncertainty.h new file mode 100644 index 0000000..e36457f --- /dev/null +++ b/src/VofTracking/Uncertainty.h @@ -0,0 +1,18 @@ +#pragma once + +#include +#include +#include + +#include "Grid/DomainInfo.h" +#include "Seed.h" + +namespace VofFlow { + struct UncertaintyArrays { + vtkSmartPointer sameEndPhaseRatio; + vtkSmartPointer stayedInPhaseRatio; + }; + + UncertaintyArrays makeUncertaintyArrays(const DomainInfo& domainInfo, const SeedPoints& seeds, + const SeedResult& seedResult); +} // namespace VofFlow diff --git a/src/VofTracking/vtk.module b/src/VofTracking/vtk.module new file mode 100644 index 0000000..039dda9 --- /dev/null +++ b/src/VofTracking/vtk.module @@ -0,0 +1,8 @@ +NAME + VofFlow::VofTracking +DEPENDS + VTK::CommonCore + VTK::FiltersCore + VTK::FiltersGeneral + VTK::ParallelMPI + VofFlow::VofDataUtils diff --git a/src/VofTracking/vtkVofTracking.cxx b/src/VofTracking/vtkVofTracking.cxx new file mode 100644 index 0000000..5f61220 --- /dev/null +++ b/src/VofTracking/vtkVofTracking.cxx @@ -0,0 +1,664 @@ +#include "vtkVofTracking.h" + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Advection.h" +#include "Boundary.h" +#include "ComponentExtraction.h" +#include "Constants.h" +#include "Grid/DomainInfo.h" +#include "Grid/GridTypes.h" +#include "Math/Vector.h" +#include "Misc/ListSearch.h" +#include "Misc/Profiling.h" +#include "Misc/VofData.h" +#include "Particles.h" +#include "Plic/CachedPlic.h" +#include "Seed.h" +#include "Uncertainty.h" +#include "VtkUtil/VtkUtilArray.h" +#include "VtkUtil/VtkUtilMPI.h" + +#define MEASURE_TIME +#include "TimeMeasure.h" + +#ifndef _WIN32 +#include +#endif + +using json = nlohmann::json; + +vtkStandardNewMacro(vtkVofTracking); + +vtkVofTracking::vtkVofTracking() + : UseThreePhase(0), + UseComponents(0), + UseTargetTimeStep(0), + InitTimeStep(0), + TargetTimeStep(0), + Refinement(0), + NeighborCorrection(0), + CellCorrection(0), + PLICCorrection(0), + IntegrationMethod(1), + IntegrationSubSteps(8), + Epsilon(0.0), + NumIterations(20), + GhostCells(4), + CutLabels(0), + CutFunction(nullptr), + BoundaryMethod(2), + OutputDataType(2), + OutputState(1), + OutputTimeMeasure(0), + MirrorXMin(0), + MirrorXMax(0), + MirrorYMin(0), + MirrorYMax(0), + MirrorZMin(0), + MirrorZMax(0), + timestepT0_(-1), + timestepT1_(-1), + timestepInit_(-1), + timestepTarget_(-1), + internalState_(NONE), + lastMTime_(0), + mpiController_(nullptr), + requiredNumGhostLevels_(4), + domainInfo_(nullptr), + seedInfo_(nullptr), + lastDataVelocity_(nullptr), + seedPoints_(nullptr), + particles_(nullptr), + oobParticles_(nullptr) { + ZoneScoped; + + internalMTime_.Modified(); + + this->SetNumberOfInputPorts(1); + this->SetNumberOfOutputPorts(4); + + mpiController_ = vtkMPIController::SafeDownCast(vtkMultiProcessController::GetGlobalController()); +#ifndef _WIN32 + // Log pid to allow attaching debugger to correct MPI rank. + if (mpiController_ != nullptr) { + vtkLog(INFO, "PID: " << getpid()); + } +#endif +} + +vtkVofTracking::~vtkVofTracking() = default; + +void vtkVofTracking::PrintSelf(ostream& os, vtkIndent indent) { + this->Superclass::PrintSelf(os, indent); +} + +vtkMTimeType vtkVofTracking::GetMTime() { + vtkMTimeType maxMTime = this->Superclass::GetMTime(); + + if (CutFunction != nullptr) { + vtkMTimeType funcMTime = CutFunction->GetMTime(); + if (funcMTime > maxMTime) { + maxMTime = funcMTime; + } + } + return maxMTime; +} + +void vtkVofTracking::Modified() { + internalMTime_.Modified(); + this->Superclass::Modified(); +} + +void vtkVofTracking::ModifiedNotInternally() { + this->Superclass::Modified(); +} + +int vtkVofTracking::FillInputPortInformation(int port, vtkInformation* info) { + ZoneScoped; + + // Note: The vtkRectilinearGrid output must use port 0. DO NOT CHANGE THIS! + // We have a vtkRectilinearGrid input for this filter and ParaView seems to automatically propagate some extent + // information if input and output types are the same. I assume, there is a bug in VTK/ParaView where maybe port 0 + // is hardcoded for some if this automatic information propagation or something similar. So, this does not work, + // when using another port than 0. Even worse, it seems to have no effect when trying to manually set the correct + // output extent information manually on another output port. The result of this is that ParaView will call + // RequestData multiple times with a different FROM_OUTPUT_PORT() request value, when running the filter in + // parallel using MPI. + if (port == 0) { + info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkRectilinearGrid"); + return 1; + } + return vtkDataSetAlgorithm::FillInputPortInformation(port, info); +} + +int vtkVofTracking::FillOutputPortInformation(int port, vtkInformation* info) { + ZoneScoped; + + if (port == 0) { + // ParaView does not call this method after parameters are changed. In addition, ParaView seems to enforce the + // data type set here, no matter if it is overwritten in RequestDataObject. Therefore, set to generic + // `vtkDataObject` instead of setting to `vtkRectilinearGrid` or `vtkImageData` depending on OutputDataType. + info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkDataObject"); + return 1; + } else if (port <= 3) { + info->Set(vtkDataObject::DATA_TYPE_NAME(), "vtkPolyData"); + return 1; + } + return vtkDataSetAlgorithm::FillOutputPortInformation(port, info); +} + +int vtkVofTracking::RequestInformation(vtkInformation* vtkNotUsed(request), vtkInformationVector** inputVector, + vtkInformationVector* vtkNotUsed(outputVector)) { + ZoneScoped; + + if (isRank0()) { + vtkLog(INFO, "RequestInformation"); + } + + vtkInformation* inInfoGrid = inputVector[0]->GetInformationObject(0); + + if (!inInfoGrid->Has(vtkStreamingDemandDrivenPipeline::TIME_STEPS())) { + vtkErrorMacro(<< "Input information has no TIME_STEPS set"); + return 0; + } + + auto numberOfTimeSteps = inInfoGrid->Length(vtkStreamingDemandDrivenPipeline::TIME_STEPS()); + if (numberOfTimeSteps <= 1) { + vtkWarningMacro(<< "Not enough input time steps for topology computation"); + } + + timeStepValues_.resize(numberOfTimeSteps); + inInfoGrid->Get(vtkStreamingDemandDrivenPipeline::TIME_STEPS(), timeStepValues_.data()); + + // Validate order of time steps, should be provided sorted from ParaView + for (std::size_t i = 1; i < timeStepValues_.size(); i++) { + if (timeStepValues_[i - 1] > timeStepValues_[i]) { + vtkErrorMacro(<< "Time steps are not in increasing order!"); + return 0; + } + } + + return 1; +} + +int vtkVofTracking::RequestUpdateExtent(vtkInformation* vtkNotUsed(request), vtkInformationVector** inputVector, + vtkInformationVector* outputVector) { + ZoneScoped; + + if (isRank0()) { + vtkLog(INFO, "RequestUpdateExtent"); + } + + vtkInformation* inInfoGrid = inputVector[0]->GetInformationObject(0); + + // set ghost level + if (GhostCells < requiredNumGhostLevels_) { + vtkWarningMacro("At least 4 ghost cells are required!"); + } + int outNumGhostLevels = hasMPI() ? std::max(requiredNumGhostLevels_, GhostCells) : 0; + for (int i = 0; i < outputVector->GetNumberOfInformationObjects(); i++) { + vtkInformation* outInfo = outputVector->GetInformationObject(i); + outNumGhostLevels = std::max(outNumGhostLevels, + outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS())); + } + inInfoGrid->Set(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(), outNumGhostLevels); + + // init and target time step + timestepInit_ = std::clamp(InitTimeStep, 0, static_cast(timeStepValues_.size() - 1)); + if (UseTargetTimeStep) { + timestepTarget_ = TargetTimeStep; + } else { + // Use UPDATE_TIME_STEP from input. Somehow it seems, that the VTK pipeline is meant to request this from the + // output port (not 100% sure about this). However, we have multiple output ports here, and it seems ParaView + // only updates one output port with the currently selected timestep from the UI, all others just receive the + // previous time step and hidden output ports get never updated. So we are required to determine the correct + // output port (could maybe be the request property FROM_OUTPUT_PORT, but don't know for sure). + // Instead, it seems that ParaView copies the value from the correct output port to the input port before + // RequestUpdateExtent runs. Therefore, we are just reading from the input, not sure if this is meant to be + // done this way, but it works so far. Side note: This only works in RequestUpdateExtent as expected. Later in + // RequestData only the time step we request here from the streaming pipeline is returned. + double targetTime = inInfoGrid->Get(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP()); + timestepTarget_ = static_cast(VofFlow::findClosestIndex(targetTime, timeStepValues_)); + } + if (timestepTarget_ < timestepInit_) { + timestepTarget_ = timestepInit_; + vtkWarningMacro("TargetTimeStep is smaller than InitTimeStep! Changed TargetTimeStep to InitTimeStep."); + } + timestepTarget_ = std::clamp(timestepTarget_, 0, static_cast(timeStepValues_.size() - 1)); + + // State Setup + // Restart algorithm if internal MTime changed or current timestep is above target. If we have not yet reached the + // target time, or only the target time was increased, we can continue with advection. Otherwise, output is reached + // or an update not changing the internal MTime happened (i.e. the cut function changed). + if (internalMTime_.GetMTime() != lastMTime_ || timestepT1_ > timestepTarget_ || internalState_ == NONE) { + lastMTime_ = internalMTime_.GetMTime(); + internalState_ = INIT; + timestepT0_ = timestepInit_; + timestepT1_ = timestepInit_; + } else if (timestepT1_ < timestepTarget_) { + internalState_ = ADVECTION; + timestepT0_ = timestepT1_; + timestepT1_++; + } else { + internalState_ = OUTPUT; + timestepT0_ = timestepT1_; + } + + // Set time step for RequestData call + if (timestepT1_ <= timestepTarget_) { + inInfoGrid->Set(vtkStreamingDemandDrivenPipeline::UPDATE_TIME_STEP(), timeStepValues_[timestepT1_]); + } + + return 1; +} + +int vtkVofTracking::RequestDataObject(vtkInformation* vtkNotUsed(request), + vtkInformationVector** vtkNotUsed(inputVector), vtkInformationVector* outputVector) { + ZoneScoped; + + if (isRank0()) { + vtkLog(INFO, "RequestDataObject"); + } + + for (int i = 0; i < this->GetNumberOfOutputPorts(); ++i) { + vtkInformation* info = outputVector->GetInformationObject(i); + vtkDataSet* output = vtkDataSet::SafeDownCast(info->Get(vtkDataObject::DATA_OBJECT())); + if (i == 0) { + if (OutputDataType == 1) { + if (!output || !output->IsA("vtkRectilinearGrid")) { + output = vtkRectilinearGrid::New(); + info->Set(vtkDataObject::DATA_OBJECT(), output); + output->Delete(); + } + } else if (OutputDataType == 2) { + if (!output || !output->IsA("vtkImageData")) { + output = vtkImageData::New(); + info->Set(vtkDataObject::DATA_OBJECT(), output); + output->Delete(); + } + } else { + vtkErrorMacro(<< "Invalid OutputDataType!"); + return 0; + } + } else { + if (!output) { + output = vtkPolyData::New(); + info->Set(vtkDataObject::DATA_OBJECT(), output); + output->Delete(); + } + } + } + return 1; +} + +int vtkVofTracking::RequestData(vtkInformation* request, vtkInformationVector** inputVector, + vtkInformationVector* outputVector) { + ZoneScoped; + + if (isRank0()) { + auto state = "RequestData " + toString(internalState_) + " T0: " + std::to_string(timestepT0_) + + " T1: " + std::to_string(timestepT1_); + ZoneName(state.c_str(), state.size()); + vtkLog(INFO, << state); + } + + vtkInformation* inInfoGrid = inputVector[0]->GetInformationObject(0); + + vtkRectilinearGrid* inputGrid = vtkRectilinearGrid::SafeDownCast(inInfoGrid->Get(vtkDataObject::DATA_OBJECT())); + if (!inputGrid) { + vtkErrorMacro(<< "Input grid missing!"); + return 0; + } + + // Mapping in plugin xml: 0 = f; 1 = f3 (optional); 2 = norm (optional); 3 = velocity + VofFlow::VofData inDataVof{ + GetInputArrayToProcess(UseThreePhase ? 1 : 0, inputGrid), + UseThreePhase ? GetInputArrayToProcess(0, inputGrid) : nullptr, + UseThreePhase ? GetInputArrayToProcess(2, inputGrid) : nullptr, + }; + vtkSmartPointer inDataVelocity = GetInputArrayToProcess(3, inputGrid); + if (inDataVof.vof1st == nullptr || + (UseThreePhase && (inDataVof.vof2nd == nullptr || inDataVof.vof2ndNorm == nullptr))) { + vtkErrorMacro(<< "VoF input array missing!"); + return 0; + } + + vtkSmartPointer inDataComponents1 = nullptr; + vtkSmartPointer inDataComponents2 = nullptr; + if (UseComponents) { + inDataComponents1 = GetInputArrayToProcess(UseThreePhase ? 5 : 4, inputGrid); + inDataComponents2 = UseThreePhase ? GetInputArrayToProcess(4, inputGrid) : nullptr; + if (inDataComponents1 == nullptr || (UseThreePhase && inDataComponents2 == nullptr)) { + vtkErrorMacro(<< "Component input array missing!"); + return 0; + } + } + + vtkInformation* outInfo0 = outputVector->GetInformationObject(0); + vtkInformation* outInfo1 = outputVector->GetInformationObject(1); + vtkInformation* outInfo2 = outputVector->GetInformationObject(2); + vtkInformation* outInfo3 = outputVector->GetInformationObject(3); + vtkDataSet* outputGrid = vtkDataSet::SafeDownCast(outInfo0->Get(vtkDataObject::DATA_OBJECT())); + vtkPolyData* outputSeeds = vtkPolyData::SafeDownCast(outInfo1->Get(vtkDataObject::DATA_OBJECT())); + vtkPolyData* outputParticles = vtkPolyData::SafeDownCast(outInfo2->Get(vtkDataObject::DATA_OBJECT())); + vtkPolyData* outputBoundary = vtkPolyData::SafeDownCast(outInfo3->Get(vtkDataObject::DATA_OBJECT())); + + try { + // === Init ============================================================ + if (internalState_ == INIT) { + TIME_MEASURE_RESET(mpiController_); + TIME_MEASURE_SYNC_START(TimeMeasure::Name::All); + + TIME_MEASURE_SYNC_START(TimeMeasure::Name::Seeding); + + domainInfo_ = std::make_unique(inputGrid, mpiController_); + seedInfo_ = std::make_unique(*domainInfo_, Refinement); + + if (OutputDataType == 2 && !domainInfo_->isUniform()) { + vtkErrorMacro("Input grid is not uniform. Cannot output vtkImageData!"); + } + + seedPoints_ = generateSeedPoints(*domainInfo_, *seedInfo_, inDataVof, Epsilon, NumIterations); + + auto numSeeds = seedPoints_->points->GetNumberOfTuples(); + if (domainInfo_->isParallel()) { + numSeeds = VofFlow::ReduceSum(mpiController_, numSeeds, 0); + } + if (isRank0()) { + vtkLog(INFO, "Number of seeded points: " << numSeeds); + } + + int processId = domainInfo_->isParallel() ? mpiController_->GetLocalProcessId() : -1; + particles_ = VofFlow::initParticles(*seedPoints_, processId); + oobParticles_ = std::make_unique(); + + TIME_MEASURE_END(TimeMeasure::Name::Seeding); + } else { + // Validate domain info by checking if extent did not change over time. + VofFlow::extent_t inExtent{}; + inputGrid->GetExtent(inExtent.data()); + if (domainInfo_->localExtent() != inExtent) { + vtkErrorMacro(<< "Input data extent changed! This is not supported!"); + return 0; + } + } + + // Globally shared PLIC cache + VofFlow::CachedPlic plicCache(*domainInfo_, inDataVof, Epsilon, NumIterations); + + // === Advection ============================================================ + if (internalState_ == ADVECTION) { + TIME_MEASURE_SYNC_START(TimeMeasure::Name::Advection); + + if (inDataVelocity != nullptr) { + TIME_MEASURE_SYNC_START(TimeMeasure::Name::Advection_Init); + double dt = timeStepValues_[timestepT1_] - timeStepValues_[timestepT0_]; + const auto numTimeSteps = timestepT1_ - timestepInit_; + VofFlow::mirrorInfo_t mirrorInfo{MirrorXMin != 0, MirrorXMax != 0, MirrorYMin != 0, MirrorYMax != 0, + MirrorZMin != 0, MirrorZMax != 0}; + + std::vector particlesPosOld = particles_->position; + + TIME_MEASURE_END(TimeMeasure::Name::Advection_Init); + + TIME_MEASURE_SYNC_START(TimeMeasure::Name::Advection_Advection); + advectParticles(*domainInfo_, lastDataVelocity_, inDataVelocity, particles_->position, dt, + IntegrationMethod, IntegrationSubSteps, mirrorInfo); + TIME_MEASURE_END(TimeMeasure::Name::Advection_Advection); + + TIME_MEASURE_SYNC_START(TimeMeasure::Name::Advection_OOB); + VofFlow::extractOutOfBoundsParticles(*particles_, particlesPosOld, *oobParticles_, *domainInfo_, + numTimeSteps); + TIME_MEASURE_END(TimeMeasure::Name::Advection_OOB); + + TIME_MEASURE_SYNC_START(TimeMeasure::Name::Advection_Corrector); + VofFlow::correctParticles(*domainInfo_, *particles_, particlesPosOld, inDataVof, NeighborCorrection, + CellCorrection, PLICCorrection, Epsilon, plicCache); + particlesPosOld.clear(); + TIME_MEASURE_END(TimeMeasure::Name::Advection_Corrector); + + // Check if particles switched phase. If all three correctors are enabled, there is no phase change. + if (NeighborCorrection == 0 || CellCorrection == 0 || PLICCorrection == 0) { + TIME_MEASURE_SYNC_START(TimeMeasure::Name::Advection_PhaseChange); + VofFlow::checkPhaseChanged(*particles_, plicCache, numTimeSteps); + TIME_MEASURE_END(TimeMeasure::Name::Advection_PhaseChange); + } + + if (domainInfo_->isParallel()) { + TIME_MEASURE_SYNC_START(TimeMeasure::Name::Advection_Exchange); + VofFlow::exchangeParticles(*particles_, *domainInfo_, mpiController_); + TIME_MEASURE_END(TimeMeasure::Name::Advection_Exchange); + } + } + + TIME_MEASURE_END(TimeMeasure::Name::Advection); + } + + // === Output ============================================================ + if (internalState_ == OUTPUT) { + + TIME_MEASURE_SYNC_START(TimeMeasure::Name::Components); + + // === Component Extraction ============================================================ + VofFlow::ComponentResult components1; + VofFlow::ComponentResult components2; + VofFlow::ComponentExtractor extractor(*domainInfo_, mpiController_); + components1 = UseComponents ? extractor.checkComponents(inDataComponents1) + : extractor.extractComponents(inDataVof.vof1st); + if (inDataVof.vof2nd != nullptr) { + components2 = UseComponents ? extractor.checkComponents(inDataComponents2) + : extractor.extractComponents(inDataVof.vof2nd); + VofFlow::appendArrayName(components1.gridLabels, " [1]"); + VofFlow::appendArrayName(components2.gridLabels, " [2]"); + } + + TIME_MEASURE_END(TimeMeasure::Name::Components); + + TIME_MEASURE_SYNC_START(TimeMeasure::Name::ParticleLabeling); + + // === Particle Labels ============================================================ + std::vector particleLabels = VofFlow::labelAdvectedParticles(*domainInfo_, components1, components2, + *particles_, CutLabels ? CutFunction : nullptr); + + std::vector particleTargetPhases = VofFlow::getParticlePhase(plicCache, *particles_); + + TIME_MEASURE_END(TimeMeasure::Name::ParticleLabeling); + + TIME_MEASURE_SYNC_START(TimeMeasure::Name::SeedLabeling); + + // === Seed Labels ============================================================ + const auto seedResult = VofFlow::transferParticleDataToSeeds(*seedPoints_, *particles_, particleLabels, + particleTargetPhases, *oobParticles_, domainInfo_->isParallel() ? mpiController_ : nullptr); + + TIME_MEASURE_END(TimeMeasure::Name::SeedLabeling); + + TIME_MEASURE_SYNC_START(TimeMeasure::Name::Boundary); + + // === Boundaries ============================================================ + VofFlow::BoundarySeedPoints boundarySeedPoints(seedPoints_->seedIdx, seedResult.labels); + VofFlow::BoundarySeedPoints neighborBoundarySeedPoints; + if (domainInfo_->isParallel()) { + neighborBoundarySeedPoints = + exchangeBoundarySeeds(*domainInfo_, *seedInfo_, boundarySeedPoints, mpiController_); + } + + vtkSmartPointer boundaries = VofFlow::generateBoundary(*domainInfo_, *seedInfo_, + boundarySeedPoints, neighborBoundarySeedPoints, BoundaryMethod, false); + + TIME_MEASURE_END(TimeMeasure::Name::Boundary); + + TIME_MEASURE_SYNC_START(TimeMeasure::Name::Output); + + // === Output Data ============================================================ + vtkSmartPointer particles = particles_->toPolyData(); + particles->GetPointData()->AddArray( + VofFlow::createVtkArray(VofFlow::ArrayNames::LABELS, particleLabels)); + + float uncertainty = std::reduce(particles_->uncertainty.cbegin(), particles_->uncertainty.cend()); + if (domainInfo_->isParallel()) { + uncertainty = VofFlow::ReduceSum(mpiController_, uncertainty, 0); + } + if (isRank0()) { + vtkLog(INFO, "Uncertainty: " << uncertainty); + } + + vtkSmartPointer glyphFilter = vtkSmartPointer::New(); + + glyphFilter->SetInputData(VofFlow::makePolyData(*seedPoints_, seedResult)); + glyphFilter->Update(); + outputSeeds->ShallowCopy(glyphFilter->GetOutput()); + + glyphFilter->SetInputData(particles); + glyphFilter->Update(); + outputParticles->ShallowCopy(glyphFilter->GetOutput()); + + outputBoundary->ShallowCopy(boundaries); + + if (OutputDataType == 1) { + vtkRectilinearGrid::SafeDownCast(outputGrid)->CopyStructure(domainInfo_->grid()); + } else if (OutputDataType == 2) { + auto* img = vtkImageData::SafeDownCast(outputGrid); + const auto& ext = domainInfo_->localExtent(); + const auto& dims = domainInfo_->cellDims(); + const auto& bounds = domainInfo_->localBounds(); + img->SetExtent(ext[0], ext[1], ext[2], ext[3], ext[4], ext[5]); + img->SetOrigin(bounds[0], bounds[2], bounds[4]); + img->SetSpacing((bounds[1] - bounds[0]) / dims[0], (bounds[3] - bounds[2]) / dims[1], + (bounds[5] - bounds[4]) / dims[2]); + } else { + vtkErrorMacro(<< "Invalid OutputDataType!"); + return 0; + } + + outputGrid->GetCellData()->AddArray(components1.gridLabels); + if (inDataVof.vof2nd != nullptr) { + outputGrid->GetCellData()->AddArray(components2.gridLabels); + } + + const auto uncertaintyArrays = VofFlow::makeUncertaintyArrays(*domainInfo_, *seedPoints_, seedResult); + outputGrid->GetCellData()->AddArray(uncertaintyArrays.sameEndPhaseRatio); + outputGrid->GetCellData()->AddArray(uncertaintyArrays.stayedInPhaseRatio); + + // Crop ghost cells from output + outputGrid->Crop(domainInfo_->localZeroExtent().data()); + + if (OutputState) { + auto state = VofFlow::createVtkStringArrayFromString("State", stateJson()); + outputGrid->GetFieldData()->AddArray(state); + outputSeeds->GetFieldData()->AddArray(state); + outputParticles->GetFieldData()->AddArray(state); + outputBoundary->GetFieldData()->AddArray(state); + } + + TIME_MEASURE_END(TimeMeasure::Name::Output); + + TIME_MEASURE_SYNC_END(TimeMeasure::Name::All); + if (OutputTimeMeasure) { + json j; + TIME_MEASURE_JSON(j); + auto timeMeasure = VofFlow::createVtkStringArrayFromString("TimeMeasure", j.dump()); + outputGrid->GetFieldData()->AddArray(timeMeasure); + outputSeeds->GetFieldData()->AddArray(timeMeasure); + outputParticles->GetFieldData()->AddArray(timeMeasure); + outputBoundary->GetFieldData()->AddArray(timeMeasure); + } + } + } catch (const std::exception& ex) { + vtkErrorMacro(<< ex.what()); + return 0; + } + + // Store time step velocity for integration + lastDataVelocity_ = inDataVelocity; + + // Set pipeline state + if (internalState_ == OUTPUT) { + request->Remove(vtkStreamingDemandDrivenPipeline::CONTINUE_EXECUTING()); + } else { + request->Set(vtkStreamingDemandDrivenPipeline::CONTINUE_EXECUTING(), 1); + } + return 1; +} + +std::string vtkVofTracking::stateJson() const { + json parameters; + parameters["UseThreePhase"] = UseThreePhase; + parameters["UseComponents"] = UseComponents; + parameters["UseTargetTimeStep"] = UseTargetTimeStep; + parameters["InitTimeStep"] = InitTimeStep; + parameters["TargetTimeStep"] = TargetTimeStep; + parameters["Refinement"] = Refinement; + parameters["NeighborCorrection"] = NeighborCorrection; + parameters["CellCorrection"] = CellCorrection; + parameters["PLICCorrection"] = PLICCorrection; + parameters["IntegrationMethod"] = IntegrationMethod; + parameters["IntegrationSubSteps"] = IntegrationSubSteps; + parameters["Epsilon"] = Epsilon; + parameters["NumIterations"] = NumIterations; + parameters["GhostCells"] = GhostCells; + parameters["CutLabels"] = CutLabels; + parameters["CutFunction"] = json(); // TODO + parameters["BoundaryMethod"] = BoundaryMethod; + parameters["OutputDataType"] = OutputDataType; + parameters["OutputState"] = OutputState; + parameters["OutputTimeMeasure"] = OutputTimeMeasure; + parameters["MirrorXMin"] = MirrorXMin; + parameters["MirrorXMax"] = MirrorXMax; + parameters["MirrorYMin"] = MirrorYMin; + parameters["MirrorYMax"] = MirrorYMax; + parameters["MirrorZMin"] = MirrorZMin; + parameters["MirrorZMax"] = MirrorZMax; + json domain; + domain["LocalExtent"] = domainInfo_->localExtent(); + domain["LocalZeroExtent"] = domainInfo_->localZeroExtent(); + domain["GlobalExtent"] = domainInfo_->globalExtent(); + domain["LocalBounds"] = domainInfo_->localBounds(); + domain["LocalZeroBounds"] = domainInfo_->localZeroBounds(); + domain["GlobalBounds"] = domainInfo_->globalBounds(); + domain["CellDims"] = domainInfo_->cellDims(); + domain["IsUniform"] = domainInfo_->isUniform(); + domain["IsParallel"] = domainInfo_->isParallel(); + json mpi; + if (domainInfo_->isParallel()) { + mpi["NumberOfProcesses"] = mpiController_->GetNumberOfProcesses(); + mpi["LocalProcessId"] = mpiController_->GetLocalProcessId(); + } + json time; + time["TimeStepInit"] = timestepInit_; + time["TimeStepTarget"] = timestepTarget_; + + json state; + state["parameters"] = parameters; + state["domain"] = domain; + state["mpi"] = mpi; + state["time"] = time; + + return state.dump(); +} diff --git a/src/VofTracking/vtkVofTracking.h b/src/VofTracking/vtkVofTracking.h new file mode 100644 index 0000000..4673adf --- /dev/null +++ b/src/VofTracking/vtkVofTracking.h @@ -0,0 +1,231 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include "VofTrackingModule.h" + +class vtkDataArray; +class vtkImplicitFunction; + +namespace VofFlow { + class DomainInfo; + class SeedCoordInfo; + struct SeedPoints; + struct Particles; + struct OutOfBoundsParticles; +} // namespace VofFlow + +#define vtkSetMacroNoInternalModified(name, type) \ + virtual void Set##name(type _arg) { \ + if (this->name != _arg) { \ + this->name = _arg; \ + this->ModifiedNotInternally(); \ + } \ + } + +class VOFTRACKING_EXPORT vtkVofTracking : public vtkDataSetAlgorithm { +public: + static vtkVofTracking* New(); + vtkTypeMacro(vtkVofTracking, vtkDataSetAlgorithm); + void PrintSelf(ostream& os, vtkIndent indent) override; + + // ParaView GUI + vtkGetMacro(UseThreePhase, int); + vtkSetMacro(UseThreePhase, int); + + vtkGetMacro(UseComponents, int); + vtkSetMacro(UseComponents, int); + + vtkGetMacro(UseTargetTimeStep, int); + vtkSetMacro(UseTargetTimeStep, int); + + vtkGetMacro(InitTimeStep, int); + vtkSetMacro(InitTimeStep, int); + + vtkGetMacro(TargetTimeStep, int); + vtkSetMacro(TargetTimeStep, int); + + vtkGetMacro(Refinement, int); + vtkSetMacro(Refinement, int); + + vtkGetMacro(NeighborCorrection, int); + vtkSetMacro(NeighborCorrection, int); + + vtkGetMacro(CellCorrection, int); + vtkSetMacro(CellCorrection, int); + + vtkGetMacro(PLICCorrection, int); + vtkSetMacro(PLICCorrection, int); + + vtkGetMacro(IntegrationMethod, int); + vtkSetMacro(IntegrationMethod, int); + + vtkGetMacro(IntegrationSubSteps, int); + vtkSetMacro(IntegrationSubSteps, int); + + vtkGetMacro(Epsilon, double); + vtkSetMacro(Epsilon, double); + + vtkGetMacro(NumIterations, int); + vtkSetMacro(NumIterations, int); + + vtkGetMacro(GhostCells, int); + vtkSetMacro(GhostCells, int); + + vtkGetMacro(CutLabels, int); + vtkSetMacroNoInternalModified(CutLabels, int); + + vtkGetMacro(CutFunction, vtkImplicitFunction*); + vtkSetMacroNoInternalModified(CutFunction, vtkImplicitFunction*); + + vtkGetMacro(BoundaryMethod, int); + vtkSetMacro(BoundaryMethod, int); + + vtkGetMacro(OutputDataType, int); + vtkSetMacro(OutputDataType, int); + + vtkGetMacro(OutputState, int); + vtkSetMacro(OutputState, int); + + vtkGetMacro(OutputTimeMeasure, int); + vtkSetMacro(OutputTimeMeasure, int); + + vtkGetMacro(MirrorXMin, int); + vtkSetMacro(MirrorXMin, int); + + vtkGetMacro(MirrorXMax, int); + vtkSetMacro(MirrorXMax, int); + + vtkGetMacro(MirrorYMin, int); + vtkSetMacro(MirrorYMin, int); + + vtkGetMacro(MirrorYMax, int); + vtkSetMacro(MirrorYMax, int); + + vtkGetMacro(MirrorZMin, int); + vtkSetMacro(MirrorZMin, int); + + vtkGetMacro(MirrorZMax, int); + vtkSetMacro(MirrorZMax, int); + + // Override GetMTime to also check for CutFunction changes. + vtkMTimeType GetMTime() override; + + // We want to be able to change some of our filter properties without needing to rerun our full algorithm. + // Therefore, we use internalMTime_ to track state in addition the MTime of this vtkObject. To be notified of all + // external changes we need to override Modified. In addition, ModifiedNotInternally only updates the MTime of the + // vtkObject to let ParaView know there was a change, but without updating the internal MTime to not rerun the + // algorithm. + void Modified() override; + void ModifiedNotInternally(); + + vtkVofTracking(const vtkVofTracking&) = delete; + void operator=(const vtkVofTracking&) = delete; + +protected: + vtkVofTracking(); + ~vtkVofTracking() override; + + int FillInputPortInformation(int port, vtkInformation* info) override; + int FillOutputPortInformation(int port, vtkInformation* info) override; + int RequestInformation(vtkInformation* request, vtkInformationVector** inputVector, + vtkInformationVector* outputVector) override; + int RequestUpdateExtent(vtkInformation* request, vtkInformationVector** inputVector, + vtkInformationVector* outputVector) override; + int RequestDataObject(vtkInformation* request, vtkInformationVector** inputVector, + vtkInformationVector* outputVector) override; + int RequestData(vtkInformation* request, vtkInformationVector** inputVector, + vtkInformationVector* outputVector) override; + +private: + [[nodiscard]] inline bool hasMPI() const { + return mpiController_ != nullptr && mpiController_->GetCommunicator() != nullptr; + } + + [[nodiscard]] inline bool isRank0() const { + return mpiController_ == nullptr || mpiController_->GetLocalProcessId() == 0; + } + + [[nodiscard]] std::string stateJson() const; + + enum State { + NONE = 0, + INIT, + ADVECTION, + OUTPUT, + }; + + static inline std::string toString(State s) { + switch (s) { + case NONE: + return "NONE"; + case INIT: + return "INIT"; + case ADVECTION: + return "ADVECTION"; + case OUTPUT: + return "OUTPUT"; + } + return "INVALID STATE"; + } + + // ParaView GUI + int UseThreePhase; + int UseComponents; + int UseTargetTimeStep; + int InitTimeStep; + int TargetTimeStep; + int Refinement; + int NeighborCorrection; + int CellCorrection; + int PLICCorrection; + int IntegrationMethod; + int IntegrationSubSteps; + double Epsilon; + int NumIterations; + int GhostCells; + int CutLabels; + vtkImplicitFunction* CutFunction; + int BoundaryMethod; + int OutputDataType; + int OutputState; + int OutputTimeMeasure; + int MirrorXMin; + int MirrorXMax; + int MirrorYMin; + int MirrorYMax; + int MirrorZMin; + int MirrorZMax; + + // Time steps + state + std::vector timeStepValues_; + int timestepT0_; + int timestepT1_; + int timestepInit_; + int timestepTarget_; + State internalState_; + vtkTimeStamp internalMTime_; + vtkMTimeType lastMTime_; + + // MPI / domain context + vtkMPIController* mpiController_; + const int requiredNumGhostLevels_; + std::unique_ptr domainInfo_; + std::unique_ptr seedInfo_; + + // Data from last time step + vtkSmartPointer lastDataVelocity_; + + // Particles + std::unique_ptr seedPoints_; + std::unique_ptr particles_; + std::unique_ptr oobParticles_; +}; diff --git a/vcpkg.json b/vcpkg.json new file mode 100644 index 0000000..ad6a642 --- /dev/null +++ b/vcpkg.json @@ -0,0 +1,8 @@ +{ + "name": "vofflow", + "version": "1.0.0", + "dependencies": [ + "cgal", + "nlohmann-json" + ] +}