From 7193f8ba0d6858638faf57ae82ce2406376ca6f2 Mon Sep 17 00:00:00 2001 From: Andrea Paris Date: Fri, 27 Feb 2026 12:26:15 +0100 Subject: [PATCH 01/12] Flow analysis tool --- tsd/src/tsd/core/algorithms/vort.h | 288 +++++++++++++++ tsd/src/tsd/io/CMakeLists.txt | 1 + tsd/src/tsd/io/importers/import_VTI.cpp | 83 ++++- tsd/src/tsd/io/procedural.hpp | 1 + .../tsd/io/procedural/computeVorticity.cpp | 336 ++++++++++++++++++ .../tsd/io/procedural/computeVorticity.hpp | 37 ++ tsd/src/tsd/ui/imgui/Application.cpp | 11 + tsd/src/tsd/ui/imgui/Application.h | 2 + tsd/src/tsd/ui/imgui/CMakeLists.txt | 1 + .../tsd/ui/imgui/modals/VorticityDialog.cpp | 135 +++++++ tsd/src/tsd/ui/imgui/modals/VorticityDialog.h | 31 ++ 11 files changed, 916 insertions(+), 10 deletions(-) create mode 100644 tsd/src/tsd/core/algorithms/vort.h create mode 100644 tsd/src/tsd/io/procedural/computeVorticity.cpp create mode 100644 tsd/src/tsd/io/procedural/computeVorticity.hpp create mode 100644 tsd/src/tsd/ui/imgui/modals/VorticityDialog.cpp create mode 100644 tsd/src/tsd/ui/imgui/modals/VorticityDialog.h diff --git a/tsd/src/tsd/core/algorithms/vort.h b/tsd/src/tsd/core/algorithms/vort.h new file mode 100644 index 000000000..886192ad7 --- /dev/null +++ b/tsd/src/tsd/core/algorithms/vort.h @@ -0,0 +1,288 @@ +// Copyright 2025-2026 NVIDIA Corporation +// SPDX-License-Identifier: Apache-2.0 + +#pragma once + +#include +#include +#include +#include + +// 3x3 symmetric eigenvalue solver — analytical (Cardano/trigonometric method) +// Returns eigenvalues sorted largest→smallest in eigs[0..2] +static void eigen3sym(double A[3][3], double eigs[3]) +{ + double p1 = A[0][1] * A[0][1] + A[0][2] * A[0][2] + A[1][2] * A[1][2]; + if (p1 == 0.0) { + // Diagonal matrix — eigenvalues are the diagonal entries + eigs[0] = A[0][0]; + eigs[1] = A[1][1]; + eigs[2] = A[2][2]; + if (eigs[0] < eigs[1]) + std::swap(eigs[0], eigs[1]); + if (eigs[1] < eigs[2]) + std::swap(eigs[1], eigs[2]); + if (eigs[0] < eigs[1]) + std::swap(eigs[0], eigs[1]); + return; + } + double q = (A[0][0] + A[1][1] + A[2][2]) / 3.0; + double p2 = (A[0][0] - q) * (A[0][0] - q) + (A[1][1] - q) * (A[1][1] - q) + + (A[2][2] - q) * (A[2][2] - q) + 2 * p1; + double p = std::sqrt(p2 / 6.0); + double B[3][3]; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + B[i][j] = (A[i][j] - (i == j ? q : 0.0)) / p; + double r = (B[0][0] * (B[1][1] * B[2][2] - B[1][2] * B[2][1]) + - B[0][1] * (B[1][0] * B[2][2] - B[1][2] * B[2][0]) + + B[0][2] * (B[1][0] * B[2][1] - B[1][1] * B[2][0])) + / 2.0; + double phi; + if (r <= -1.0) + phi = M_PI / 3.0; + else if (r >= 1.0) + phi = 0.0; + else + phi = std::acos(r) / 3.0; + eigs[0] = q + 2 * p * std::cos(phi); + eigs[2] = q + 2 * p * std::cos(phi + 2.0943951023931953); // phi + 2*pi/3 + eigs[1] = 3 * q - eigs[0] - eigs[2]; +} + +// lambda2: second (middle) eigenvalue of S^2 + O^2 +// J is the velocity gradient Jacobian; uses no external dependencies +static double l2(double J[3][3]) +{ + double S[3][3], O[3][3]; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) { + S[i][j] = 0.5 * (J[i][j] + J[j][i]); + O[i][j] = S[i][j] - J[j][i]; // O = S - J^T = 0.5*(J - J^T) + } + double M[3][3] = {}; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + for (int k = 0; k < 3; k++) + M[i][j] += S[i][k] * S[k][j] + O[i][k] * O[k][j]; + double eigs[3]; + eigen3sym(M, eigs); // sorted largest→smallest; eigs[1] is the middle + return eigs[1]; +} + +// Q-criterion: 0.5 * (||O||_F^2 - ||S||_F^2) +static double q_crit(double J[3][3]) +{ + double S[3][3], O[3][3]; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) { + S[i][j] = 0.5 * (J[i][j] + J[j][i]); + O[i][j] = S[i][j] - J[j][i]; + } + double trO2 = 0, trS2 = 0; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) { + trO2 += O[i][j] * O[i][j]; + trS2 += S[i][j] * S[i][j]; + } + return 0.5 * (trO2 - trS2); +} + +// Forward declarations + +inline void gradX(std::vector &u, + std::vector &uGrad, + double dx, + size_t nx, + size_t ny, + size_t nz); +inline void gradY(std::vector &v, + std::vector &vGrad, + double dy, + size_t nx, + size_t ny, + size_t nz); +inline void gradZ(std::vector &w, + std::vector &wGrad, + std::vector &z_, + size_t nx, + size_t ny, + size_t nz); +inline void vort(std::vector &u, + std::vector &v, + std::vector &w, + std::vector &x_, + std::vector &y_, + std::vector &z_, + std::vector &vorticity, + std::vector &helicity, + std::vector &lambda2, + std::vector &qCriterion, + size_t nx, + size_t ny, + size_t nz); + +// 2nd-order one-sided BCs, 2nd-order central FD in interior, in x +inline void gradX(std::vector &u, + std::vector &uGrad, + double dx, + size_t nx, + size_t ny, + size_t nz) +{ +#pragma omp parallel for + for (auto z = 0; z < (int)nz; ++z) { + auto off = z * nx * ny; + for (auto row = 0; row < (int)ny; ++row) { + // Left boundary: 2nd-order one-sided forward + uGrad[off + row * nx + 0] = + (-3.0 * u[off + row * nx + 0] + 4.0 * u[off + row * nx + 1] + - u[off + row * nx + 2]) + / (2.0 * dx); + // Right boundary: 2nd-order one-sided backward + uGrad[off + row * nx + nx - 1] = + (3.0 * u[off + row * nx + nx - 1] - 4.0 * u[off + row * nx + nx - 2] + + u[off + row * nx + nx - 3]) + / (2.0 * dx); + // Interior: 2nd-order central + for (auto col = 1; col < (int)nx - 1; ++col) { + uGrad[off + row * nx + col] = + (u[off + row * nx + col + 1] - u[off + row * nx + col - 1]) + / (2.0 * dx); + } + } + } +} + +// 2nd-order one-sided BCs, 2nd-order central FD in interior, in y +inline void gradY(std::vector &v, + std::vector &vGrad, + double dy, + size_t nx, + size_t ny, + size_t nz) +{ +#pragma omp parallel for + for (auto z = 0; z < (int)nz; ++z) { + auto off = z * nx * ny; + for (auto col = 0; col < (int)nx; ++col) { + // Bottom boundary (row=0): 2nd-order one-sided forward + vGrad[0 * nx + col + off] = + (-3.0 * v[0 * nx + col + off] + 4.0 * v[1 * nx + col + off] + - v[2 * nx + col + off]) + / (2.0 * dy); + // Top boundary (row=ny-1): 2nd-order one-sided backward + vGrad[(ny - 1) * nx + col + off] = + (3.0 * v[(ny - 1) * nx + col + off] + - 4.0 * v[(ny - 2) * nx + col + off] + + v[(ny - 3) * nx + col + off]) + / (2.0 * dy); + // Interior: 2nd-order central + for (auto row = 1; row < (int)ny - 1; ++row) { + vGrad[row * nx + col + off] = + (v[(row + 1) * nx + col + off] - v[(row - 1) * nx + col + off]) + / (2.0 * dy); + } + } + } +} + +// Non-uniform spacing in z (one-sided FD at boundaries) +inline void gradZ(std::vector &w, + std::vector &wGrad, + std::vector &z_, + size_t nx, + size_t ny, + size_t nz) +{ + size_t off = nx * ny; +#pragma omp parallel for + for (auto row = 0; row < (int)ny; ++row) { + for (auto col = 0; col < (int)nx; ++col) { + wGrad[(nz - 1) * off + row * nx + col] = + (w[(nz - 1) * off + row * nx + col] + - w[(nz - 2) * off + row * nx + col]) + / (z_[nz - 1] - z_[nz - 2]); + wGrad[0 * off + row * nx + col] = + (w[1 * off + row * nx + col] - w[0 * off + row * nx + col]) + / (z_[1] - z_[0]); + // Interior: 2nd-order central + for (auto z = 1; z < (int)nz - 1; ++z) { + wGrad[z * off + row * nx + col] = + (w[(z + 1) * off + row * nx + col] + - w[(z - 1) * off + row * nx + col]) + / (z_[z + 1] - z_[z - 1]); + } + } + } +} + +inline void vort(std::vector &u, + std::vector &v, + std::vector &w, + std::vector &x_, + std::vector &y_, + std::vector &z_, + std::vector &vorticity, + std::vector &helicity, + std::vector &lambda2, + std::vector &qCriterion, + size_t nx, + size_t ny, + size_t nz) +{ + auto len = u.size(); + std::vector dux(len); + std::vector dvx(len); + std::vector dwx(len); + std::vector duy(len); + std::vector dvy(len); + std::vector dwy(len); + std::vector duz(len); + std::vector dvz(len); + std::vector dwz(len); + + double dx = x_.at(1) - x_.at(0); + double dy = y_.at(1) - y_.at(0); + + gradX(u, dux, dx, nx, ny, nz); + gradX(v, dvx, dx, nx, ny, nz); + gradX(w, dwx, dx, nx, ny, nz); + + gradY(u, duy, dy, nx, ny, nz); + gradY(v, dvy, dy, nx, ny, nz); + gradY(w, dwy, dy, nx, ny, nz); + + gradZ(u, duz, z_, nx, ny, nz); + gradZ(v, dvz, z_, nx, ny, nz); + gradZ(w, dwz, z_, nx, ny, nz); + + for (size_t i = 0; i < len; ++i) { + vorticity[i] = + std::sqrt((dwy[i] - dvz[i]) * (dwy[i] - dvz[i]) + + (duz[i] - dwx[i]) * (duz[i] - dwx[i]) + + (dvx[i] - duy[i]) * (dvx[i] - duy[i])); + + helicity[i] = std::abs( + ((dwy[i] - dvz[i]) * u[i]) + ((duz[i] - dwx[i]) * v[i]) + + ((dvx[i] - duy[i]) * w[i])); + + double velMag = std::sqrt(u[i] * u[i] + v[i] * v[i] + w[i] * w[i]); + if (velMag != 0 && vorticity[i] != 0) { + helicity[i] /= (2 * velMag * vorticity[i]); + } else { + helicity[i] = 0; + } + } + +#pragma omp parallel for + for (size_t i = 0; i < len; ++i) { + double J[3][3] = {{dux[i], duy[i], duz[i]}, + {dvx[i], dvy[i], dvz[i]}, + {dwx[i], dwy[i], dwz[i]}}; + + lambda2[i] = -std::min(l2(J), 0.0); + qCriterion[i] = std::max(q_crit(J), 0.0); + } + std::cout << "[vort] Vortical variables computed" << std::endl; +} diff --git a/tsd/src/tsd/io/CMakeLists.txt b/tsd/src/tsd/io/CMakeLists.txt index 15d199ed9..d839e4799 100644 --- a/tsd/src/tsd/io/CMakeLists.txt +++ b/tsd/src/tsd/io/CMakeLists.txt @@ -42,6 +42,7 @@ PRIVATE importers/import_VTP.cpp importers/import_VTU.cpp importers/import_XYZDP.cpp + procedural/computeVorticity.cpp procedural/generate_cylinders.cpp procedural/generate_default_lights.cpp procedural/generate_hdri_dome.cpp diff --git a/tsd/src/tsd/io/importers/import_VTI.cpp b/tsd/src/tsd/io/importers/import_VTI.cpp index 1270238ca..e0d9b1f89 100644 --- a/tsd/src/tsd/io/importers/import_VTI.cpp +++ b/tsd/src/tsd/io/importers/import_VTI.cpp @@ -51,25 +51,88 @@ SpatialFieldRef import_VTI(Scene &scene, const char *filepath) field->setParameter("origin", float3(origin[0], origin[1], origin[2])); field->setParameter("spacing", float3(spacing[0], spacing[1], spacing[2])); - // --- Write point data array --- + // --- Write point data arrays --- vtkPointData *pointData = grid->GetPointData(); for (uint32_t i = 0; i < pointData->GetNumberOfArrays(); ++i) { vtkDataArray *array = pointData->GetArray(i); int numComponents = array->GetNumberOfComponents(); - if (numComponents > 1) { + + if (numComponents == 1) { + auto a = makeArray3DFromVTK( + scene, array, dims[0], dims[1], dims[2], "[import_VTI]"); + field->setParameterObject("data", *a); + break; + } else if (numComponents == 3) { + // Split into 3 scalar SpatialFields: {name}_x, {name}_y, {name}_z + const char *arrName = array->GetName(); + std::string baseName = (arrName && arrName[0] != '\0') + ? std::string(arrName) + : fileOf(filepath); + + std::string nameX = baseName + "_x"; + std::string nameY = baseName + "_y"; + std::string nameZ = baseName + "_z"; + + size_t n = + (size_t)dims[0] * (size_t)dims[1] * (size_t)dims[2]; + + auto arrX = scene.createArray(ANARI_FLOAT32, dims[0], dims[1], dims[2]); + auto arrY = scene.createArray(ANARI_FLOAT32, dims[0], dims[1], dims[2]); + auto arrZ = scene.createArray(ANARI_FLOAT32, dims[0], dims[1], dims[2]); + + float *pX = arrX->mapAs(); + float *pY = arrY->mapAs(); + float *pZ = arrZ->mapAs(); + + for (vtkIdType idx = 0; idx < (vtkIdType)n; ++idx) { + double tuple[3] = {0.0, 0.0, 0.0}; + array->GetTuple(idx, tuple); + pX[idx] = (float)tuple[0]; + pY[idx] = (float)tuple[1]; + pZ[idx] = (float)tuple[2]; + } + + arrX->unmap(); + arrY->unmap(); + arrZ->unmap(); + + // Main field (_x component) + field->setName(nameX.c_str()); + field->setParameterObject("data", *arrX); + + // _y component field in the scene object pool + auto fieldY = scene.createObject( + tokens::spatial_field::structuredRegular); + fieldY->setName(nameY.c_str()); + fieldY->setParameter("origin", float3(origin[0], origin[1], origin[2])); + fieldY->setParameter( + "spacing", float3(spacing[0], spacing[1], spacing[2])); + fieldY->setParameterObject("data", *arrY); + + // _z component field in the scene object pool + auto fieldZ = scene.createObject( + tokens::spatial_field::structuredRegular); + fieldZ->setName(nameZ.c_str()); + fieldZ->setParameter("origin", float3(origin[0], origin[1], origin[2])); + fieldZ->setParameter( + "spacing", float3(spacing[0], spacing[1], spacing[2])); + fieldZ->setParameterObject("data", *arrZ); + + logStatus( + "[import_VTI] split 3-component array '%s' into '%s', '%s', '%s'", + baseName.c_str(), + nameX.c_str(), + nameY.c_str(), + nameZ.c_str()); + break; + } else { logWarning( - "[import_VTI] only single-component arrays are supported, " - "array '%s' has %d components -- only using first component", + "[import_VTI] array '%s' has %d components (only 1 or 3 are " + "supported) -- skipping", array->GetName(), numComponents); - continue; } - - auto a = makeArray3DFromVTK( - scene, array, dims[0], dims[1], dims[2], "[import_VTI]"); - field->setParameterObject("data", *a); - break; } return field; diff --git a/tsd/src/tsd/io/procedural.hpp b/tsd/src/tsd/io/procedural.hpp index d67522b9f..f92ac751e 100644 --- a/tsd/src/tsd/io/procedural.hpp +++ b/tsd/src/tsd/io/procedural.hpp @@ -4,6 +4,7 @@ #pragma once #include "tsd/core/scene/Scene.hpp" +#include "tsd/io/procedural/computeVorticity.hpp" namespace tsd::io { diff --git a/tsd/src/tsd/io/procedural/computeVorticity.cpp b/tsd/src/tsd/io/procedural/computeVorticity.cpp new file mode 100644 index 000000000..7fc59e88a --- /dev/null +++ b/tsd/src/tsd/io/procedural/computeVorticity.cpp @@ -0,0 +1,336 @@ +// Copyright 2025-2026 NVIDIA Corporation +// SPDX-License-Identifier: Apache-2.0 + +#include "computeVorticity.hpp" +// tsd_core +#include "tsd/core/ColorMapUtil.hpp" +#include "tsd/core/Logging.hpp" +#include "tsd/core/algorithms/vort.h" +// nanovdb +#include +// anari +#include +// std +#include +#include + +namespace tsd::io { + +using namespace tsd::core; + +// --------------------------------------------------------------------------- +// Internal helpers +// --------------------------------------------------------------------------- + +struct FieldData +{ + std::vector values; + std::vector x, y, z; + size_t nx{0}, ny{0}, nz{0}; +}; + +static bool extractStructuredRegular( + Scene &scene, SpatialField *field, FieldData &out) +{ + auto *p = field->parameter("data"); + if (!p || !anari::isArray(p->value().type())) { + logError("[computeVorticity] structuredRegular field has no 'data' array"); + return false; + } + + auto arr = scene.getObject(p->value().getAsObjectIndex()); + if (!arr) { + logError("[computeVorticity] structuredRegular 'data' array is invalid"); + return false; + } + + out.nx = arr->dim(0); + out.ny = arr->dim(1); + out.nz = arr->dim(2); + + if (out.nx < 2 || out.ny < 2 || out.nz < 2) { + logError("[computeVorticity] field dimensions must be >= 2 in each axis"); + return false; + } + + if (arr->elementType() != ANARI_FLOAT32) { + logError( + "[computeVorticity] structuredRegular field element type must be " + "ANARI_FLOAT32 (got %d)", + (int)arr->elementType()); + return false; + } + + auto *p_orig = field->parameter("origin"); + auto *p_spc = field->parameter("spacing"); + if (!p_orig || !p_spc) { + logError("[computeVorticity] structuredRegular field missing origin/spacing"); + return false; + } + + auto origin = p_orig->value().get(); + auto spacing = p_spc->value().get(); + + // Build coordinate arrays + out.x.resize(out.nx); + out.y.resize(out.ny); + out.z.resize(out.nz); + for (size_t i = 0; i < out.nx; ++i) + out.x[i] = origin.x + i * (double)spacing.x; + for (size_t j = 0; j < out.ny; ++j) + out.y[j] = origin.y + j * (double)spacing.y; + for (size_t k = 0; k < out.nz; ++k) + out.z[k] = origin.z + k * (double)spacing.z; + + // Convert float → double + size_t total = out.nx * out.ny * out.nz; + out.values.resize(total); + const float *src = arr->dataAs(); + for (size_t i = 0; i < total; ++i) + out.values[i] = src[i]; + + return true; +} + +static bool extractNanoVDB(Scene &scene, SpatialField *field, FieldData &out) +{ + auto *p = field->parameter("data"); + if (!p || !anari::isArray(p->value().type())) { + logError("[computeVorticity] nanovdb field has no 'data' array"); + return false; + } + + auto arr = scene.getObject(p->value().getAsObjectIndex()); + if (!arr || !arr->data()) { + logError("[computeVorticity] nanovdb 'data' array is invalid"); + return false; + } + + const uint8_t *rawData = static_cast(arr->data()); + const auto *meta = + reinterpret_cast(rawData); + + if (meta->gridType() != nanovdb::GridType::Float) { + logError( + "[computeVorticity] nanovdb field must be float32 (GridType::Float)"); + return false; + } + + const auto *grid = + reinterpret_cast *>(rawData); + auto bbox = grid->indexBBox(); + auto lo = bbox.min(); + auto hi = bbox.max(); + + out.nx = (size_t)(hi[0] - lo[0] + 1); + out.ny = (size_t)(hi[1] - lo[1] + 1); + out.nz = (size_t)(hi[2] - lo[2] + 1); + + if (out.nx < 2 || out.ny < 2 || out.nz < 2) { + logError("[computeVorticity] nanovdb field dimensions must be >= 2"); + return false; + } + + // Derive world-space coordinates from the grid map + const auto &map = grid->map(); + auto worldLo = map.applyMap(nanovdb::Vec3d(lo[0], lo[1], lo[2])); + auto worldHi = map.applyMap(nanovdb::Vec3d(hi[0], hi[1], hi[2])); + + double dx = (out.nx > 1) ? (worldHi[0] - worldLo[0]) / (out.nx - 1) : 1.0; + double dy = (out.ny > 1) ? (worldHi[1] - worldLo[1]) / (out.ny - 1) : 1.0; + double dz = (out.nz > 1) ? (worldHi[2] - worldLo[2]) / (out.nz - 1) : 1.0; + + out.x.resize(out.nx); + out.y.resize(out.ny); + out.z.resize(out.nz); + for (size_t i = 0; i < out.nx; ++i) + out.x[i] = worldLo[0] + i * dx; + for (size_t j = 0; j < out.ny; ++j) + out.y[j] = worldLo[1] + j * dy; + for (size_t k = 0; k < out.nz; ++k) + out.z[k] = worldLo[2] + k * dz; + + // Rasterize into a dense buffer using the accessor + size_t total = out.nx * out.ny * out.nz; + out.values.resize(total, 0.0); + auto acc = grid->getAccessor(); + + for (int k = lo[2]; k <= hi[2]; ++k) { + size_t kk = (size_t)(k - lo[2]); + for (int j = lo[1]; j <= hi[1]; ++j) { + size_t jj = (size_t)(j - lo[1]); + for (int i = lo[0]; i <= hi[0]; ++i) { + size_t ii = (size_t)(i - lo[0]); + out.values[kk * out.ny * out.nx + jj * out.nx + ii] = + (double)acc.getValue(nanovdb::Coord(i, j, k)); + } + } + } + + return true; +} + +static bool extractFieldData(Scene &scene, SpatialField *field, FieldData &out) +{ + if (!field) { + logError("[computeVorticity] null SpatialField pointer"); + return false; + } + + if (field->subtype() == tokens::spatial_field::structuredRegular) { + return extractStructuredRegular(scene, field, out); + } else if (field->subtype() == tokens::spatial_field::nanovdb) { + return extractNanoVDB(scene, field, out); + } else { + logError( + "[computeVorticity] unsupported SpatialField subtype '%s'; " + "only structuredRegular and nanovdb are supported", + field->subtype().c_str()); + return false; + } +} + +static VolumeRef makeOutputVolume(Scene &scene, + const std::string &name, + const std::vector &data, + size_t nx, + size_t ny, + size_t nz, + const math::float3 &origin, + const math::float3 &spacing, + LayerNodeRef location) +{ + // Create output SpatialField + auto field = + scene.createObject(tokens::spatial_field::structuredRegular); + field->setName(name.c_str()); + field->setParameter("origin", origin); + field->setParameter("spacing", spacing); + + // Fill output array + auto dataArr = scene.createArray(ANARI_FLOAT32, nx, ny, nz); + float *dst = dataArr->mapAs(); + for (size_t i = 0; i < data.size(); ++i) + dst[i] = (float)data[i]; + dataArr->unmap(); + field->setParameterObject("data", *dataArr); + + // Compute value range + float2 valueRange = field->computeValueRange(); + + // Create Volume node in the scene tree + auto tx = scene.insertChildTransformNode( + location ? location : scene.defaultLayer()->root()); + + auto [inst, vol] = scene.insertNewChildObjectNode( + tx, tokens::volume::transferFunction1D); + vol->setName(name.c_str()); + vol->setParameterObject("value", *field); + vol->setParameter("valueRange", ANARI_FLOAT32_BOX1, &valueRange); + + // Default colormap + auto colorArr = scene.createArray(ANARI_FLOAT32_VEC4, 256); + colorArr->setData(makeDefaultColorMap(256).data()); + vol->setParameterObject("color", *colorArr); + + return vol; +} + +// --------------------------------------------------------------------------- +// Public API +// --------------------------------------------------------------------------- + +VorticityResult computeVorticity(Scene &scene, + SpatialField *u, + SpatialField *v, + SpatialField *w, + LayerNodeRef location, + VorticityOptions opts) +{ + VorticityResult result; + + logStatus("[computeVorticity] extracting field data..."); + + FieldData uData, vData, wData; + if (!extractFieldData(scene, u, uData)) + return result; + if (!extractFieldData(scene, v, vData)) + return result; + if (!extractFieldData(scene, w, wData)) + return result; + + // Verify matching dimensions + if (uData.nx != vData.nx || uData.nx != wData.nx || uData.ny != vData.ny + || uData.ny != wData.ny || uData.nz != vData.nz + || uData.nz != wData.nz) { + logError( + "[computeVorticity] U/V/W fields have mismatched dimensions: " + "U=(%zu,%zu,%zu) V=(%zu,%zu,%zu) W=(%zu,%zu,%zu)", + uData.nx, + uData.ny, + uData.nz, + vData.nx, + vData.ny, + vData.nz, + wData.nx, + wData.ny, + wData.nz); + return result; + } + + size_t nx = uData.nx, ny = uData.ny, nz = uData.nz; + size_t total = nx * ny * nz; + + logStatus( + "[computeVorticity] computing vortical quantities on %zux%zux%zu grid...", + nx, + ny, + nz); + + // Allocate output buffers + std::vector vorticityBuf(total, 0.0); + std::vector helicityBuf(total, 0.0); + std::vector lambda2Buf(total, 0.0); + std::vector qCritBuf(total, 0.0); + + // Run the vorticity computation + vort(uData.values, + vData.values, + wData.values, + uData.x, + uData.y, + uData.z, + vorticityBuf, + helicityBuf, + lambda2Buf, + qCritBuf, + nx, + ny, + nz); + + logStatus("[computeVorticity] creating output volumes..."); + + // Derive origin and spacing from the U field coordinates + math::float3 origin{(float)uData.x[0], (float)uData.y[0], (float)uData.z[0]}; + math::float3 spacing{(float)(uData.x[1] - uData.x[0]), + (float)(uData.y[1] - uData.y[0]), + (float)(uData.z[1] - uData.z[0])}; + + if (opts.lambda2) + result.lambda2 = makeOutputVolume( + scene, "lambda2", lambda2Buf, nx, ny, nz, origin, spacing, location); + if (opts.qCriterion) + result.qCriterion = makeOutputVolume( + scene, "q_criterion", qCritBuf, nx, ny, nz, origin, spacing, location); + if (opts.vorticity) + result.vorticity = makeOutputVolume( + scene, "vorticity", vorticityBuf, nx, ny, nz, origin, spacing, location); + if (opts.helicity) + result.helicity = makeOutputVolume( + scene, "helicity", helicityBuf, nx, ny, nz, origin, spacing, location); + + logStatus("[computeVorticity] done."); + return result; +} + +} // namespace tsd::io diff --git a/tsd/src/tsd/io/procedural/computeVorticity.hpp b/tsd/src/tsd/io/procedural/computeVorticity.hpp new file mode 100644 index 000000000..65bce5688 --- /dev/null +++ b/tsd/src/tsd/io/procedural/computeVorticity.hpp @@ -0,0 +1,37 @@ +// Copyright 2025-2026 NVIDIA Corporation +// SPDX-License-Identifier: Apache-2.0 + +#pragma once + +#include "tsd/core/scene/Scene.hpp" +#include "tsd/core/scene/objects/SpatialField.hpp" +#include "tsd/core/scene/objects/Volume.hpp" + +namespace tsd::io { + +using namespace tsd::core; + +struct VorticityOptions +{ + bool lambda2{true}; + bool qCriterion{true}; + bool vorticity{true}; + bool helicity{false}; +}; + +struct VorticityResult +{ + VolumeRef lambda2; + VolumeRef qCriterion; + VolumeRef vorticity; + VolumeRef helicity; +}; + +VorticityResult computeVorticity(Scene &scene, + SpatialField *u, + SpatialField *v, + SpatialField *w, + LayerNodeRef location = {}, + VorticityOptions opts = {}); + +} // namespace tsd::io diff --git a/tsd/src/tsd/ui/imgui/Application.cpp b/tsd/src/tsd/ui/imgui/Application.cpp index 943273d65..bb5809783 100644 --- a/tsd/src/tsd/ui/imgui/Application.cpp +++ b/tsd/src/tsd/ui/imgui/Application.cpp @@ -128,6 +128,7 @@ anari_viewer::WindowArray Application::setupWindows() m_offlineRenderModal = std::make_unique(this); m_fileDialog = std::make_unique(this); m_exportNanoVDBFileDialog = std::make_unique(this); + m_vorticityDialog = std::make_unique(this); m_applicationName = SDL_GetWindowTitle(sdlWindow()); updateWindowTitle(); @@ -187,6 +188,11 @@ void Application::uiFrameStart() modalActive = true; } + if (m_vorticityDialog->visible()) { + m_vorticityDialog->renderUI(); + modalActive = true; + } + if (ImGui::IsKeyChordPressed(ImGuiMod_Ctrl | ImGuiMod_Shift | ImGuiKey_S)) this->getFilenameFromDialog(m_filenameToSaveNextFrame, true); else if (ImGui::IsKeyChordPressed(ImGuiMod_Ctrl | ImGuiMod_Alt | ImGuiKey_S)) @@ -335,6 +341,11 @@ void Application::uiMainMenuBar_Tools() ImGui::Separator(); + if (ImGui::MenuItem("Flow Analysis...")) + m_vorticityDialog->show(); + + ImGui::Separator(); + if (ImGui::BeginMenu("TSD Device")) { if (tsdDeviceIsSetup()) { if (ImGui::MenuItem("Disable")) diff --git a/tsd/src/tsd/ui/imgui/Application.h b/tsd/src/tsd/ui/imgui/Application.h index 3170337af..3dbfab9a8 100644 --- a/tsd/src/tsd/ui/imgui/Application.h +++ b/tsd/src/tsd/ui/imgui/Application.h @@ -10,6 +10,7 @@ #include "modals/ExportNanoVDBFileDialog.h" #include "modals/ImportFileDialog.h" #include "modals/OfflineRenderModal.h" +#include "modals/VorticityDialog.h" // tsd_app #include "tsd/app/Core.h" // tsd_core @@ -117,6 +118,7 @@ class Application : public anari_viewer::Application std::unique_ptr m_offlineRenderModal; std::unique_ptr m_fileDialog; std::unique_ptr m_exportNanoVDBFileDialog; + std::unique_ptr m_vorticityDialog; tsd::core::DataTree m_settings; diff --git a/tsd/src/tsd/ui/imgui/CMakeLists.txt b/tsd/src/tsd/ui/imgui/CMakeLists.txt index 13125c5cd..dda23bb21 100644 --- a/tsd/src/tsd/ui/imgui/CMakeLists.txt +++ b/tsd/src/tsd/ui/imgui/CMakeLists.txt @@ -14,6 +14,7 @@ project_sources(PRIVATE modals/OfflineRenderModal.cpp modals/ImportFileDialog.cpp modals/Modal.cpp + modals/VorticityDialog.cpp windows/Animations.cpp windows/CameraPoses.cpp windows/DatabaseEditor.cpp diff --git a/tsd/src/tsd/ui/imgui/modals/VorticityDialog.cpp b/tsd/src/tsd/ui/imgui/modals/VorticityDialog.cpp new file mode 100644 index 000000000..8116e9be9 --- /dev/null +++ b/tsd/src/tsd/ui/imgui/modals/VorticityDialog.cpp @@ -0,0 +1,135 @@ +// Copyright 2025-2026 NVIDIA Corporation +// SPDX-License-Identifier: Apache-2.0 + +#include "VorticityDialog.h" +// tsd_io +#include "tsd/io/procedural/computeVorticity.hpp" +// tsd_core +#include "tsd/core/ObjectPool.hpp" +// tsd_ui_imgui +#include "tsd/ui/imgui/Application.h" +#include "tsd/ui/imgui/modals/BlockingTaskModal.h" +// imgui +#include "imgui.h" +// std +#include +#include + +namespace tsd::ui::imgui { + +VorticityDialog::VorticityDialog(Application *app) + : Modal(app, "Flow Analysis") +{} + +void VorticityDialog::buildUI() +{ + using namespace tsd::core; + + auto *core = appCore(); + auto &scene = core->tsd.scene; + + // Collect all SpatialFields in the scene + std::vector fields; + std::vector fieldNames; + + foreach_item_const(scene.objectDB().field, [&](const SpatialField *sf) { + if (sf) { + fields.push_back(const_cast(sf)); + fieldNames.push_back(sf->name()); + } + }); + + // Helper to build the combo label list + std::vector labels; + labels.reserve(fields.size()); + for (auto &n : fieldNames) + labels.push_back(n.c_str()); + + const int nFields = (int)fields.size(); + + ImGui::Text("Select velocity component fields:"); + ImGui::Spacing(); + + // Clamp indices in case fields changed + if (m_uIdx >= nFields) + m_uIdx = -1; + if (m_vIdx >= nFields) + m_vIdx = -1; + if (m_wIdx >= nFields) + m_wIdx = -1; + + auto comboField = [&](const char *label, int &idx) { + const char *preview = (idx >= 0 && idx < nFields) ? labels[idx] : "(none)"; + if (ImGui::BeginCombo(label, preview)) { + for (int i = 0; i < nFields; ++i) { + bool selected = (i == idx); + if (ImGui::Selectable(labels[i], selected)) + idx = i; + if (selected) + ImGui::SetItemDefaultFocus(); + } + ImGui::EndCombo(); + } + }; + + comboField("U (x-velocity)", m_uIdx); + comboField("V (y-velocity)", m_vIdx); + comboField("W (z-velocity)", m_wIdx); + + ImGui::Spacing(); + ImGui::Separator(); + ImGui::Text("Outputs to compute:"); + + ImGui::Checkbox("Lambda2", &m_computeLambda2); + ImGui::SameLine(); + ImGui::Checkbox("Q-criterion", &m_computeQ); + ImGui::SameLine(); + ImGui::Checkbox("Vorticity magnitude", &m_computeVorticity); + ImGui::SameLine(); + ImGui::Checkbox("Helicity", &m_computeHelicity); + + ImGui::Spacing(); + ImGui::Separator(); + + if (ImGui::Button("Cancel") || ImGui::IsKeyDown(ImGuiKey_Escape)) + this->hide(); + + ImGui::SameLine(); + + bool canRun = (m_uIdx >= 0 && m_vIdx >= 0 && m_wIdx >= 0 + && m_uIdx != m_vIdx && m_uIdx != m_wIdx && m_vIdx != m_wIdx); + + ImGui::BeginDisabled(!canRun); + + if (ImGui::Button("Compute")) { + this->hide(); + + SpatialField *uField = fields[m_uIdx]; + SpatialField *vField = fields[m_vIdx]; + SpatialField *wField = fields[m_wIdx]; + + tsd::io::VorticityOptions opts; + opts.lambda2 = m_computeLambda2; + opts.qCriterion = m_computeQ; + opts.vorticity = m_computeVorticity; + opts.helicity = m_computeHelicity; + + auto doCompute = [&scene, uField, vField, wField, opts]() { + tsd::io::computeVorticity(scene, uField, vField, wField, {}, opts); + }; + + m_app->showTaskModal(doCompute, "Please Wait: Computing Flow Analysis..."); + } + + ImGui::EndDisabled(); + + if (!canRun && nFields > 0) { + ImGui::SameLine(); + ImGui::TextDisabled("(select 3 distinct fields)"); + } else if (nFields == 0) { + ImGui::SameLine(); + ImGui::TextDisabled("(no SpatialFields in scene)"); + } +} + +} // namespace tsd::ui::imgui diff --git a/tsd/src/tsd/ui/imgui/modals/VorticityDialog.h b/tsd/src/tsd/ui/imgui/modals/VorticityDialog.h new file mode 100644 index 000000000..9351b734a --- /dev/null +++ b/tsd/src/tsd/ui/imgui/modals/VorticityDialog.h @@ -0,0 +1,31 @@ +// Copyright 2025-2026 NVIDIA Corporation +// SPDX-License-Identifier: Apache-2.0 + +#pragma once + +#include "Modal.h" +// std +#include +#include + +namespace tsd::ui::imgui { + +struct VorticityDialog : public Modal +{ + VorticityDialog(Application *app); + ~VorticityDialog() override = default; + + void buildUI() override; + + private: + int m_uIdx{-1}; // index into m_fields list + int m_vIdx{-1}; + int m_wIdx{-1}; + + bool m_computeLambda2{true}; + bool m_computeQ{true}; + bool m_computeVorticity{true}; + bool m_computeHelicity{false}; +}; + +} // namespace tsd::ui::imgui From 143a783f52b8ec6237809ba55695e6825da68d6a Mon Sep 17 00:00:00 2001 From: Andrea Paris Date: Fri, 27 Feb 2026 12:36:33 +0100 Subject: [PATCH 02/12] avoid data copy --- tsd/src/tsd/core/algorithms/vort.h | 207 ++++++++++-------- .../tsd/io/procedural/computeVorticity.cpp | 170 +++++++------- 2 files changed, 212 insertions(+), 165 deletions(-) diff --git a/tsd/src/tsd/core/algorithms/vort.h b/tsd/src/tsd/core/algorithms/vort.h index 886192ad7..b81c99a89 100644 --- a/tsd/src/tsd/core/algorithms/vort.h +++ b/tsd/src/tsd/core/algorithms/vort.h @@ -88,43 +88,52 @@ static double q_crit(double J[3][3]) return 0.5 * (trO2 - trS2); } -// Forward declarations +// --------------------------------------------------------------------------- +// Gradient helpers — float* input, double* output (scratch). +// Each gradient function reads float velocity values and writes double-precision +// gradient values so that finite differences retain full precision. +// +// Boundary scheme: 2nd-order one-sided stencil. +// Interior scheme: 2nd-order central difference. +// --------------------------------------------------------------------------- -inline void gradX(std::vector &u, - std::vector &uGrad, +// Forward declarations +inline void gradX(const float *u, + double *uGrad, double dx, size_t nx, size_t ny, size_t nz); -inline void gradY(std::vector &v, - std::vector &vGrad, +inline void gradY(const float *v, + double *vGrad, double dy, size_t nx, size_t ny, size_t nz); -inline void gradZ(std::vector &w, - std::vector &wGrad, - std::vector &z_, +inline void gradZ(const float *w, + double *wGrad, + const double *z_, size_t nx, size_t ny, size_t nz); -inline void vort(std::vector &u, - std::vector &v, - std::vector &w, - std::vector &x_, - std::vector &y_, - std::vector &z_, - std::vector &vorticity, - std::vector &helicity, - std::vector &lambda2, - std::vector &qCriterion, +inline void vort(const float *u, + const float *v, + const float *w, + const double *x_, + const double *y_, + const double *z_, + float *vorticity, + float *helicity, + float *lambda2, + float *qCriterion, size_t nx, size_t ny, size_t nz); -// 2nd-order one-sided BCs, 2nd-order central FD in interior, in x -inline void gradX(std::vector &u, - std::vector &uGrad, +// --------------------------------------------------------------------------- + +inline void gradX(const float *u, + double *uGrad, double dx, size_t nx, size_t ny, @@ -132,7 +141,7 @@ inline void gradX(std::vector &u, { #pragma omp parallel for for (auto z = 0; z < (int)nz; ++z) { - auto off = z * nx * ny; + auto off = (size_t)z * nx * ny; for (auto row = 0; row < (int)ny; ++row) { // Left boundary: 2nd-order one-sided forward uGrad[off + row * nx + 0] = @@ -147,16 +156,16 @@ inline void gradX(std::vector &u, // Interior: 2nd-order central for (auto col = 1; col < (int)nx - 1; ++col) { uGrad[off + row * nx + col] = - (u[off + row * nx + col + 1] - u[off + row * nx + col - 1]) + ((double)u[off + row * nx + col + 1] + - (double)u[off + row * nx + col - 1]) / (2.0 * dx); } } } } -// 2nd-order one-sided BCs, 2nd-order central FD in interior, in y -inline void gradY(std::vector &v, - std::vector &vGrad, +inline void gradY(const float *v, + double *vGrad, double dy, size_t nx, size_t ny, @@ -164,7 +173,7 @@ inline void gradY(std::vector &v, { #pragma omp parallel for for (auto z = 0; z < (int)nz; ++z) { - auto off = z * nx * ny; + auto off = (size_t)z * nx * ny; for (auto col = 0; col < (int)nx; ++col) { // Bottom boundary (row=0): 2nd-order one-sided forward vGrad[0 * nx + col + off] = @@ -180,17 +189,19 @@ inline void gradY(std::vector &v, // Interior: 2nd-order central for (auto row = 1; row < (int)ny - 1; ++row) { vGrad[row * nx + col + off] = - (v[(row + 1) * nx + col + off] - v[(row - 1) * nx + col + off]) + ((double)v[(row + 1) * nx + col + off] + - (double)v[(row - 1) * nx + col + off]) / (2.0 * dy); } } } } -// Non-uniform spacing in z (one-sided FD at boundaries) -inline void gradZ(std::vector &w, - std::vector &wGrad, - std::vector &z_, +// Non-uniform z spacing: 1st-order one-sided at boundaries, 2nd-order central +// in the interior. +inline void gradZ(const float *w, + double *wGrad, + const double *z_, size_t nx, size_t ny, size_t nz) @@ -199,90 +210,110 @@ inline void gradZ(std::vector &w, #pragma omp parallel for for (auto row = 0; row < (int)ny; ++row) { for (auto col = 0; col < (int)nx; ++col) { + // Top boundary wGrad[(nz - 1) * off + row * nx + col] = - (w[(nz - 1) * off + row * nx + col] - - w[(nz - 2) * off + row * nx + col]) + ((double)w[(nz - 1) * off + row * nx + col] + - (double)w[(nz - 2) * off + row * nx + col]) / (z_[nz - 1] - z_[nz - 2]); + // Bottom boundary wGrad[0 * off + row * nx + col] = - (w[1 * off + row * nx + col] - w[0 * off + row * nx + col]) + ((double)w[1 * off + row * nx + col] + - (double)w[0 * off + row * nx + col]) / (z_[1] - z_[0]); // Interior: 2nd-order central for (auto z = 1; z < (int)nz - 1; ++z) { wGrad[z * off + row * nx + col] = - (w[(z + 1) * off + row * nx + col] - - w[(z - 1) * off + row * nx + col]) + ((double)w[(z + 1) * off + row * nx + col] + - (double)w[(z - 1) * off + row * nx + col]) / (z_[z + 1] - z_[z - 1]); } } } } -inline void vort(std::vector &u, - std::vector &v, - std::vector &w, - std::vector &x_, - std::vector &y_, - std::vector &z_, - std::vector &vorticity, - std::vector &helicity, - std::vector &lambda2, - std::vector &qCriterion, +// --------------------------------------------------------------------------- +// vort — compute vortical quantities from float velocity fields. +// +// Inputs: u, v, w — float velocity components (x,y,z), row-major [z][y][x] +// x_, y_, z_ — double coordinate arrays (length nx, ny, nz) +// Outputs: vorticity, helicity, lambda2, qCriterion — float*, written in-place. +// Any output pointer may be null; null outputs are simply skipped. +// If both lambda2 and qCriterion are null the eigenvalue loop is skipped. +// +// Internal gradient arrays are double-precision scratch space (9 × N doubles). +// --------------------------------------------------------------------------- +inline void vort(const float *u, + const float *v, + const float *w, + const double *x_, + const double *y_, + const double *z_, + float *vorticity, + float *helicity, + float *lambda2, + float *qCriterion, size_t nx, size_t ny, size_t nz) { - auto len = u.size(); - std::vector dux(len); - std::vector dvx(len); - std::vector dwx(len); - std::vector duy(len); - std::vector dvy(len); - std::vector dwy(len); - std::vector duz(len); - std::vector dvz(len); - std::vector dwz(len); + const size_t len = nx * ny * nz; + + // 9 double-precision gradient scratch arrays (shared across u/v/w) + std::vector dux(len), dvx(len), dwx(len); + std::vector duy(len), dvy(len), dwy(len); + std::vector duz(len), dvz(len), dwz(len); - double dx = x_.at(1) - x_.at(0); - double dy = y_.at(1) - y_.at(0); + const double dx = x_[1] - x_[0]; + const double dy = y_[1] - y_[0]; - gradX(u, dux, dx, nx, ny, nz); - gradX(v, dvx, dx, nx, ny, nz); - gradX(w, dwx, dx, nx, ny, nz); + gradX(u, dux.data(), dx, nx, ny, nz); + gradX(v, dvx.data(), dx, nx, ny, nz); + gradX(w, dwx.data(), dx, nx, ny, nz); - gradY(u, duy, dy, nx, ny, nz); - gradY(v, dvy, dy, nx, ny, nz); - gradY(w, dwy, dy, nx, ny, nz); + gradY(u, duy.data(), dy, nx, ny, nz); + gradY(v, dvy.data(), dy, nx, ny, nz); + gradY(w, dwy.data(), dy, nx, ny, nz); - gradZ(u, duz, z_, nx, ny, nz); - gradZ(v, dvz, z_, nx, ny, nz); - gradZ(w, dwz, z_, nx, ny, nz); + gradZ(u, duz.data(), z_, nx, ny, nz); + gradZ(v, dvz.data(), z_, nx, ny, nz); + gradZ(w, dwz.data(), z_, nx, ny, nz); - for (size_t i = 0; i < len; ++i) { - vorticity[i] = - std::sqrt((dwy[i] - dvz[i]) * (dwy[i] - dvz[i]) - + (duz[i] - dwx[i]) * (duz[i] - dwx[i]) - + (dvx[i] - duy[i]) * (dvx[i] - duy[i])); + // Vorticity magnitude and helicity (cheap loop, not worth parallelising over + // the eigen solver — keep separate so the OMP loop below is clean) + if (vorticity || helicity) { + for (size_t i = 0; i < len; ++i) { + const double omx = dwy[i] - dvz[i]; + const double omy = duz[i] - dwx[i]; + const double omz = dvx[i] - duy[i]; + const double omag = std::sqrt(omx * omx + omy * omy + omz * omz); - helicity[i] = std::abs( - ((dwy[i] - dvz[i]) * u[i]) + ((duz[i] - dwx[i]) * v[i]) - + ((dvx[i] - duy[i]) * w[i])); + if (vorticity) + vorticity[i] = (float)omag; - double velMag = std::sqrt(u[i] * u[i] + v[i] * v[i] + w[i] * w[i]); - if (velMag != 0 && vorticity[i] != 0) { - helicity[i] /= (2 * velMag * vorticity[i]); - } else { - helicity[i] = 0; + if (helicity) { + const double ui = u[i], vi = v[i], wi = w[i]; + const double h = std::abs(omx * ui + omy * vi + omz * wi); + const double vmag = std::sqrt(ui * ui + vi * vi + wi * wi); + helicity[i] = (vmag > 0.0 && omag > 0.0) + ? (float)(h / (2.0 * vmag * omag)) + : 0.0f; + } } } + // Lambda2 and Q-criterion (expensive: one eigensolve per voxel) + if (lambda2 || qCriterion) { #pragma omp parallel for - for (size_t i = 0; i < len; ++i) { - double J[3][3] = {{dux[i], duy[i], duz[i]}, - {dvx[i], dvy[i], dvz[i]}, - {dwx[i], dwy[i], dwz[i]}}; - - lambda2[i] = -std::min(l2(J), 0.0); - qCriterion[i] = std::max(q_crit(J), 0.0); + for (size_t i = 0; i < len; ++i) { + double J[3][3] = {{dux[i], duy[i], duz[i]}, + {dvx[i], dvy[i], dvz[i]}, + {dwx[i], dwy[i], dwz[i]}}; + if (lambda2) + lambda2[i] = (float)(-std::min(l2(J), 0.0)); + if (qCriterion) + qCriterion[i] = (float)std::max(q_crit(J), 0.0); + } } + std::cout << "[vort] Vortical variables computed" << std::endl; } diff --git a/tsd/src/tsd/io/procedural/computeVorticity.cpp b/tsd/src/tsd/io/procedural/computeVorticity.cpp index 7fc59e88a..f7de5f516 100644 --- a/tsd/src/tsd/io/procedural/computeVorticity.cpp +++ b/tsd/src/tsd/io/procedural/computeVorticity.cpp @@ -19,16 +19,26 @@ namespace tsd::io { using namespace tsd::core; // --------------------------------------------------------------------------- -// Internal helpers +// FieldData — holds a pointer to float velocity data plus coordinate arrays. +// +// For structuredRegular fields the data pointer is non-owning (points directly +// into the ANARI Array storage — no copy). For NanoVDB fields the sparse +// grid is rasterized into ownedData (float, not double) and ptr is set to +// point at that buffer. // --------------------------------------------------------------------------- struct FieldData { - std::vector values; - std::vector x, y, z; + const float *ptr{nullptr}; // non-owning (structuredRegular) or alias into ownedData + std::vector ownedData; // owns data for NanoVDB rasterization + std::vector x, y, z; // world-space coordinate arrays size_t nx{0}, ny{0}, nz{0}; }; +// --------------------------------------------------------------------------- +// Extraction helpers +// --------------------------------------------------------------------------- + static bool extractStructuredRegular( Scene &scene, SpatialField *field, FieldData &out) { @@ -71,7 +81,9 @@ static bool extractStructuredRegular( auto origin = p_orig->value().get(); auto spacing = p_spc->value().get(); - // Build coordinate arrays + // Point directly at the ANARI array's storage — no copy needed. + out.ptr = arr->dataAs(); + out.x.resize(out.nx); out.y.resize(out.ny); out.z.resize(out.nz); @@ -82,13 +94,6 @@ static bool extractStructuredRegular( for (size_t k = 0; k < out.nz; ++k) out.z[k] = origin.z + k * (double)spacing.z; - // Convert float → double - size_t total = out.nx * out.ny * out.nz; - out.values.resize(total); - const float *src = arr->dataAs(); - for (size_t i = 0; i < total; ++i) - out.values[i] = src[i]; - return true; } @@ -131,7 +136,6 @@ static bool extractNanoVDB(Scene &scene, SpatialField *field, FieldData &out) return false; } - // Derive world-space coordinates from the grid map const auto &map = grid->map(); auto worldLo = map.applyMap(nanovdb::Vec3d(lo[0], lo[1], lo[2])); auto worldHi = map.applyMap(nanovdb::Vec3d(hi[0], hi[1], hi[2])); @@ -150,9 +154,9 @@ static bool extractNanoVDB(Scene &scene, SpatialField *field, FieldData &out) for (size_t k = 0; k < out.nz; ++k) out.z[k] = worldLo[2] + k * dz; - // Rasterize into a dense buffer using the accessor - size_t total = out.nx * out.ny * out.nz; - out.values.resize(total, 0.0); + // Rasterize sparse grid into a dense float buffer (no double conversion). + out.ownedData.resize(out.nx * out.ny * out.nz, 0.0f); + out.ptr = out.ownedData.data(); auto acc = grid->getAccessor(); for (int k = lo[2]; k <= hi[2]; ++k) { @@ -161,8 +165,8 @@ static bool extractNanoVDB(Scene &scene, SpatialField *field, FieldData &out) size_t jj = (size_t)(j - lo[1]); for (int i = lo[0]; i <= hi[0]; ++i) { size_t ii = (size_t)(i - lo[0]); - out.values[kk * out.ny * out.nx + jj * out.nx + ii] = - (double)acc.getValue(nanovdb::Coord(i, j, k)); + out.ownedData[kk * out.ny * out.nx + jj * out.nx + ii] = + acc.getValue(nanovdb::Coord(i, j, k)); } } } @@ -190,9 +194,14 @@ static bool extractFieldData(Scene &scene, SpatialField *field, FieldData &out) } } -static VolumeRef makeOutputVolume(Scene &scene, +// --------------------------------------------------------------------------- +// wrapAsVolume — create a SpatialField + transferFunction1D Volume around an +// already-filled ANARI Array. The array must already be unmapped. +// --------------------------------------------------------------------------- + +static VolumeRef wrapAsVolume(Scene &scene, const std::string &name, - const std::vector &data, + ArrayRef &dataArr, size_t nx, size_t ny, size_t nz, @@ -200,25 +209,15 @@ static VolumeRef makeOutputVolume(Scene &scene, const math::float3 &spacing, LayerNodeRef location) { - // Create output SpatialField auto field = scene.createObject(tokens::spatial_field::structuredRegular); field->setName(name.c_str()); field->setParameter("origin", origin); field->setParameter("spacing", spacing); - - // Fill output array - auto dataArr = scene.createArray(ANARI_FLOAT32, nx, ny, nz); - float *dst = dataArr->mapAs(); - for (size_t i = 0; i < data.size(); ++i) - dst[i] = (float)data[i]; - dataArr->unmap(); field->setParameterObject("data", *dataArr); - // Compute value range float2 valueRange = field->computeValueRange(); - // Create Volume node in the scene tree auto tx = scene.insertChildTransformNode( location ? location : scene.defaultLayer()->root()); @@ -228,7 +227,6 @@ static VolumeRef makeOutputVolume(Scene &scene, vol->setParameterObject("value", *field); vol->setParameter("valueRange", ANARI_FLOAT32_BOX1, &valueRange); - // Default colormap auto colorArr = scene.createArray(ANARI_FLOAT32_VEC4, 256); colorArr->setData(makeDefaultColorMap(256).data()); vol->setParameterObject("color", *colorArr); @@ -259,75 +257,93 @@ VorticityResult computeVorticity(Scene &scene, if (!extractFieldData(scene, w, wData)) return result; - // Verify matching dimensions if (uData.nx != vData.nx || uData.nx != wData.nx || uData.ny != vData.ny || uData.ny != wData.ny || uData.nz != vData.nz || uData.nz != wData.nz) { logError( "[computeVorticity] U/V/W fields have mismatched dimensions: " "U=(%zu,%zu,%zu) V=(%zu,%zu,%zu) W=(%zu,%zu,%zu)", - uData.nx, - uData.ny, - uData.nz, - vData.nx, - vData.ny, - vData.nz, - wData.nx, - wData.ny, - wData.nz); + uData.nx, uData.ny, uData.nz, + vData.nx, vData.ny, vData.nz, + wData.nx, wData.ny, wData.nz); return result; } - size_t nx = uData.nx, ny = uData.ny, nz = uData.nz; - size_t total = nx * ny * nz; + const size_t nx = uData.nx, ny = uData.ny, nz = uData.nz; logStatus( "[computeVorticity] computing vortical quantities on %zux%zux%zu grid...", - nx, - ny, - nz); + nx, ny, nz); + + // Pre-allocate one ANARI float32 array per selected output and map it so + // vort() can write directly into the final storage — no intermediate buffers. + ArrayRef lambda2Arr, qCritArr, vorticityArr, helicityArr; + float *lambda2Out = nullptr, *qCritOut = nullptr; + float *vorticityOut = nullptr, *helicityOut = nullptr; - // Allocate output buffers - std::vector vorticityBuf(total, 0.0); - std::vector helicityBuf(total, 0.0); - std::vector lambda2Buf(total, 0.0); - std::vector qCritBuf(total, 0.0); - - // Run the vorticity computation - vort(uData.values, - vData.values, - wData.values, - uData.x, - uData.y, - uData.z, - vorticityBuf, - helicityBuf, - lambda2Buf, - qCritBuf, + if (opts.lambda2) { + lambda2Arr = scene.createArray(ANARI_FLOAT32, nx, ny, nz); + lambda2Out = lambda2Arr->mapAs(); + } + if (opts.qCriterion) { + qCritArr = scene.createArray(ANARI_FLOAT32, nx, ny, nz); + qCritOut = qCritArr->mapAs(); + } + if (opts.vorticity) { + vorticityArr = scene.createArray(ANARI_FLOAT32, nx, ny, nz); + vorticityOut = vorticityArr->mapAs(); + } + if (opts.helicity) { + helicityArr = scene.createArray(ANARI_FLOAT32, nx, ny, nz); + helicityOut = helicityArr->mapAs(); + } + + // Compute — vort() reads float* velocity directly, writes float* outputs. + // Null output pointers are silently skipped inside vort(). + vort(uData.ptr, + vData.ptr, + wData.ptr, + uData.x.data(), + uData.y.data(), + uData.z.data(), + vorticityOut, + helicityOut, + lambda2Out, + qCritOut, nx, ny, nz); + // Unmap output arrays before handing them to SpatialField + if (lambda2Arr) + lambda2Arr->unmap(); + if (qCritArr) + qCritArr->unmap(); + if (vorticityArr) + vorticityArr->unmap(); + if (helicityArr) + helicityArr->unmap(); + logStatus("[computeVorticity] creating output volumes..."); - // Derive origin and spacing from the U field coordinates - math::float3 origin{(float)uData.x[0], (float)uData.y[0], (float)uData.z[0]}; - math::float3 spacing{(float)(uData.x[1] - uData.x[0]), + const math::float3 origin{ + (float)uData.x[0], (float)uData.y[0], (float)uData.z[0]}; + const math::float3 spacing{(float)(uData.x[1] - uData.x[0]), (float)(uData.y[1] - uData.y[0]), (float)(uData.z[1] - uData.z[0])}; - if (opts.lambda2) - result.lambda2 = makeOutputVolume( - scene, "lambda2", lambda2Buf, nx, ny, nz, origin, spacing, location); - if (opts.qCriterion) - result.qCriterion = makeOutputVolume( - scene, "q_criterion", qCritBuf, nx, ny, nz, origin, spacing, location); - if (opts.vorticity) - result.vorticity = makeOutputVolume( - scene, "vorticity", vorticityBuf, nx, ny, nz, origin, spacing, location); - if (opts.helicity) - result.helicity = makeOutputVolume( - scene, "helicity", helicityBuf, nx, ny, nz, origin, spacing, location); + if (lambda2Arr) + result.lambda2 = wrapAsVolume( + scene, "lambda2", lambda2Arr, nx, ny, nz, origin, spacing, location); + if (qCritArr) + result.qCriterion = wrapAsVolume( + scene, "q_criterion", qCritArr, nx, ny, nz, origin, spacing, location); + if (vorticityArr) + result.vorticity = wrapAsVolume( + scene, "vorticity", vorticityArr, nx, ny, nz, origin, spacing, location); + if (helicityArr) + result.helicity = wrapAsVolume( + scene, "helicity", helicityArr, nx, ny, nz, origin, spacing, location); logStatus("[computeVorticity] done."); return result; From 08d46eb0a3d1f8df61414fd65b814613ae9407e7 Mon Sep 17 00:00:00 2001 From: Andrea Paris Date: Fri, 27 Feb 2026 12:46:40 +0100 Subject: [PATCH 03/12] better loop indexing, removed omp remains --- tsd/src/tsd/core/algorithms/vort.h | 103 ++++++------------ .../tsd/io/procedural/computeVorticity.cpp | 12 +- 2 files changed, 39 insertions(+), 76 deletions(-) diff --git a/tsd/src/tsd/core/algorithms/vort.h b/tsd/src/tsd/core/algorithms/vort.h index b81c99a89..3927634d8 100644 --- a/tsd/src/tsd/core/algorithms/vort.h +++ b/tsd/src/tsd/core/algorithms/vort.h @@ -132,6 +132,8 @@ inline void vort(const float *u, // --------------------------------------------------------------------------- +// x is the fastest index, so every (z,y) pair is a contiguous row of nx +// floats — collapse both outer loops into a single row index. inline void gradX(const float *u, double *uGrad, double dx, @@ -139,31 +141,20 @@ inline void gradX(const float *u, size_t ny, size_t nz) { -#pragma omp parallel for - for (auto z = 0; z < (int)nz; ++z) { - auto off = (size_t)z * nx * ny; - for (auto row = 0; row < (int)ny; ++row) { - // Left boundary: 2nd-order one-sided forward - uGrad[off + row * nx + 0] = - (-3.0 * u[off + row * nx + 0] + 4.0 * u[off + row * nx + 1] - - u[off + row * nx + 2]) - / (2.0 * dx); - // Right boundary: 2nd-order one-sided backward - uGrad[off + row * nx + nx - 1] = - (3.0 * u[off + row * nx + nx - 1] - 4.0 * u[off + row * nx + nx - 2] - + u[off + row * nx + nx - 3]) - / (2.0 * dx); - // Interior: 2nd-order central - for (auto col = 1; col < (int)nx - 1; ++col) { - uGrad[off + row * nx + col] = - ((double)u[off + row * nx + col + 1] - - (double)u[off + row * nx + col - 1]) - / (2.0 * dx); - } - } + const double inv2dx = 1.0 / (2.0 * dx); + const int nrows = (int)(ny * nz); + for (int r = 0; r < nrows; ++r) { + const float *ur = u + (size_t)r * nx; + double *gr = uGrad + (size_t)r * nx; + gr[0] = (-3.0 * ur[0] + 4.0 * ur[1] - ur[2]) * inv2dx; + for (size_t c = 1; c < nx - 1; ++c) + gr[c] = ((double)ur[c + 1] - (double)ur[c - 1]) * inv2dx; + gr[nx - 1] = ( 3.0 * ur[nx-1] - 4.0 * ur[nx-2] + ur[nx-3]) * inv2dx; } } +// y has stride nx. Iterate over every (z,x) pair as the column base; +// the 1D y-pass then uses vc[y*nx] with no extra offset arithmetic. inline void gradY(const float *v, double *vGrad, double dy, @@ -171,34 +162,24 @@ inline void gradY(const float *v, size_t ny, size_t nz) { -#pragma omp parallel for - for (auto z = 0; z < (int)nz; ++z) { - auto off = (size_t)z * nx * ny; - for (auto col = 0; col < (int)nx; ++col) { - // Bottom boundary (row=0): 2nd-order one-sided forward - vGrad[0 * nx + col + off] = - (-3.0 * v[0 * nx + col + off] + 4.0 * v[1 * nx + col + off] - - v[2 * nx + col + off]) - / (2.0 * dy); - // Top boundary (row=ny-1): 2nd-order one-sided backward - vGrad[(ny - 1) * nx + col + off] = - (3.0 * v[(ny - 1) * nx + col + off] - - 4.0 * v[(ny - 2) * nx + col + off] - + v[(ny - 3) * nx + col + off]) - / (2.0 * dy); - // Interior: 2nd-order central - for (auto row = 1; row < (int)ny - 1; ++row) { - vGrad[row * nx + col + off] = - ((double)v[(row + 1) * nx + col + off] - - (double)v[(row - 1) * nx + col + off]) - / (2.0 * dy); - } + const double inv2dy = 1.0 / (2.0 * dy); + const size_t slab = nx * ny; + for (int z = 0; z < (int)nz; ++z) { + for (int x = 0; x < (int)nx; ++x) { + const float *vc = v + (size_t)z * slab + x; + double *gc = vGrad + (size_t)z * slab + x; + gc[0] = (-3.0 * vc[0] + 4.0 * vc[nx] - vc[2*nx]) * inv2dy; + for (size_t y = 1; y < ny - 1; ++y) + gc[y * nx] = ((double)vc[(y+1)*nx] - (double)vc[(y-1)*nx]) * inv2dy; + gc[(ny-1) * nx] = ( 3.0 * vc[(ny-1)*nx] - 4.0 * vc[(ny-2)*nx] + vc[(ny-3)*nx]) * inv2dy; } } } // Non-uniform z spacing: 1st-order one-sided at boundaries, 2nd-order central // in the interior. +// z has stride slab=nx*ny. Every (y,x) pair gives a z-column addressable as +// wc[k*slab] — flatten both outer loops into a single index over ny*nx. inline void gradZ(const float *w, double *wGrad, const double *z_, @@ -206,28 +187,15 @@ inline void gradZ(const float *w, size_t ny, size_t nz) { - size_t off = nx * ny; -#pragma omp parallel for - for (auto row = 0; row < (int)ny; ++row) { - for (auto col = 0; col < (int)nx; ++col) { - // Top boundary - wGrad[(nz - 1) * off + row * nx + col] = - ((double)w[(nz - 1) * off + row * nx + col] - - (double)w[(nz - 2) * off + row * nx + col]) - / (z_[nz - 1] - z_[nz - 2]); - // Bottom boundary - wGrad[0 * off + row * nx + col] = - ((double)w[1 * off + row * nx + col] - - (double)w[0 * off + row * nx + col]) - / (z_[1] - z_[0]); - // Interior: 2nd-order central - for (auto z = 1; z < (int)nz - 1; ++z) { - wGrad[z * off + row * nx + col] = - ((double)w[(z + 1) * off + row * nx + col] - - (double)w[(z - 1) * off + row * nx + col]) - / (z_[z + 1] - z_[z - 1]); - } - } + const size_t slab = nx * ny; + const int ncols = (int)slab; + for (int yx = 0; yx < ncols; ++yx) { + const float *wc = w + yx; + double *gc = wGrad + yx; + gc[0] = ((double)wc[slab] - (double)wc[0]) / (z_[1] - z_[0]); + for (size_t k = 1; k < nz - 1; ++k) + gc[k * slab] = ((double)wc[(k+1)*slab] - (double)wc[(k-1)*slab]) / (z_[k+1] - z_[k-1]); + gc[(nz-1) * slab] = ((double)wc[(nz-1)*slab] - (double)wc[(nz-2)*slab]) / (z_[nz-1] - z_[nz-2]); } } @@ -278,8 +246,6 @@ inline void vort(const float *u, gradZ(v, dvz.data(), z_, nx, ny, nz); gradZ(w, dwz.data(), z_, nx, ny, nz); - // Vorticity magnitude and helicity (cheap loop, not worth parallelising over - // the eigen solver — keep separate so the OMP loop below is clean) if (vorticity || helicity) { for (size_t i = 0; i < len; ++i) { const double omx = dwy[i] - dvz[i]; @@ -303,7 +269,6 @@ inline void vort(const float *u, // Lambda2 and Q-criterion (expensive: one eigensolve per voxel) if (lambda2 || qCriterion) { -#pragma omp parallel for for (size_t i = 0; i < len; ++i) { double J[3][3] = {{dux[i], duy[i], duz[i]}, {dvx[i], dvy[i], dvz[i]}, diff --git a/tsd/src/tsd/io/procedural/computeVorticity.cpp b/tsd/src/tsd/io/procedural/computeVorticity.cpp index f7de5f516..c14f138c9 100644 --- a/tsd/src/tsd/io/procedural/computeVorticity.cpp +++ b/tsd/src/tsd/io/procedural/computeVorticity.cpp @@ -160,14 +160,12 @@ static bool extractNanoVDB(Scene &scene, SpatialField *field, FieldData &out) auto acc = grid->getAccessor(); for (int k = lo[2]; k <= hi[2]; ++k) { - size_t kk = (size_t)(k - lo[2]); for (int j = lo[1]; j <= hi[1]; ++j) { - size_t jj = (size_t)(j - lo[1]); - for (int i = lo[0]; i <= hi[0]; ++i) { - size_t ii = (size_t)(i - lo[0]); - out.ownedData[kk * out.ny * out.nx + jj * out.nx + ii] = - acc.getValue(nanovdb::Coord(i, j, k)); - } + float *row = out.ownedData.data() + + (size_t)(k - lo[2]) * out.ny * out.nx + + (size_t)(j - lo[1]) * out.nx; + for (int i = lo[0]; i <= hi[0]; ++i) + row[i - lo[0]] = acc.getValue(nanovdb::Coord(i, j, k)); } } From c0d4c0c05ceaf482d98640bc79b39dbaac8832d6 Mon Sep 17 00:00:00 2001 From: Andrea Paris Date: Fri, 27 Feb 2026 15:00:35 +0100 Subject: [PATCH 04/12] VTU vorticity using VTK library --- tsd/src/tsd/core/algorithms/vort.h | 213 +++++--- tsd/src/tsd/io/importers/import_VTU.cpp | 91 +++- .../tsd/io/procedural/computeVorticity.cpp | 465 +++++++++++++++++- 3 files changed, 658 insertions(+), 111 deletions(-) diff --git a/tsd/src/tsd/core/algorithms/vort.h b/tsd/src/tsd/core/algorithms/vort.h index 3927634d8..bcc9bebed 100644 --- a/tsd/src/tsd/core/algorithms/vort.h +++ b/tsd/src/tsd/core/algorithms/vort.h @@ -88,25 +88,52 @@ static double q_crit(double J[3][3]) return 0.5 * (trO2 - trS2); } +// --------------------------------------------------------------------------- +// 1D finite-difference helpers — float* input, double* output, stride s. +// +// grad1D: uniform spacing; 2nd-order one-sided BCs, 2nd-order central. +// inv2h = 1/(2h). +// grad1D_nu: non-uniform spacing; 1st-order one-sided BCs, 2nd-order central. +// c[]: coordinate array (length n). +// --------------------------------------------------------------------------- + +static void grad1D( + const float *fc, double *gc, size_t n, size_t s, double inv2h) +{ + gc[0] = (-3.0*fc[0] + 4.0*fc[s] - fc[2*s]) * inv2h; + for (size_t i = 1; i < n-1; ++i) + gc[i*s] = ((double)fc[(i+1)*s] - (double)fc[(i-1)*s]) * inv2h; + gc[(n-1)*s] = ( 3.0*fc[(n-1)*s] - 4.0*fc[(n-2)*s] + fc[(n-3)*s]) * inv2h; +} + +static void grad1D_nu( + const float *fc, double *gc, size_t n, size_t s, const double *c) +{ + gc[0] = ((double)fc[s] - (double)fc[0]) / (c[1] - c[0]); + for (size_t i = 1; i < n-1; ++i) + gc[i*s] = ((double)fc[(i+1)*s] - (double)fc[(i-1)*s]) / (c[i+1] - c[i-1]); + gc[(n-1)*s] = ((double)fc[(n-1)*s] - (double)fc[(n-2)*s]) / (c[n-1] - c[n-2]); +} + // --------------------------------------------------------------------------- // Gradient helpers — float* input, double* output (scratch). -// Each gradient function reads float velocity values and writes double-precision -// gradient values so that finite differences retain full precision. +// All three take a coordinate array so uniform and rectilinear grids share the +// same interface; uniform fields pass the x_/y_/z_ arrays built from +// origin+spacing in extractStructuredRegular. // -// Boundary scheme: 2nd-order one-sided stencil. -// Interior scheme: 2nd-order central difference. +// Boundary accuracy: 1st-order one-sided. Interior: 2nd-order central. // --------------------------------------------------------------------------- // Forward declarations inline void gradX(const float *u, double *uGrad, - double dx, + const double *x_, size_t nx, size_t ny, size_t nz); inline void gradY(const float *v, double *vGrad, - double dy, + const double *y_, size_t nx, size_t ny, size_t nz); @@ -116,6 +143,23 @@ inline void gradZ(const float *w, size_t nx, size_t ny, size_t nz); +inline void vort_from_jacobians(const float *u, + const float *v, + const float *w, + const double *dux, + const double *dvx, + const double *dwx, + const double *duy, + const double *dvy, + const double *dwy, + const double *duz, + const double *dvz, + const double *dwz, + float *vorticity, + float *helicity, + float *lambda2, + float *qCriterion, + size_t len); inline void vort(const float *u, const float *v, const float *w, @@ -136,44 +180,29 @@ inline void vort(const float *u, // floats — collapse both outer loops into a single row index. inline void gradX(const float *u, double *uGrad, - double dx, + const double *x_, size_t nx, size_t ny, size_t nz) { - const double inv2dx = 1.0 / (2.0 * dx); const int nrows = (int)(ny * nz); - for (int r = 0; r < nrows; ++r) { - const float *ur = u + (size_t)r * nx; - double *gr = uGrad + (size_t)r * nx; - gr[0] = (-3.0 * ur[0] + 4.0 * ur[1] - ur[2]) * inv2dx; - for (size_t c = 1; c < nx - 1; ++c) - gr[c] = ((double)ur[c + 1] - (double)ur[c - 1]) * inv2dx; - gr[nx - 1] = ( 3.0 * ur[nx-1] - 4.0 * ur[nx-2] + ur[nx-3]) * inv2dx; - } + for (int r = 0; r < nrows; ++r) + grad1D_nu(u + (size_t)r*nx, uGrad + (size_t)r*nx, nx, 1, x_); } // y has stride nx. Iterate over every (z,x) pair as the column base; -// the 1D y-pass then uses vc[y*nx] with no extra offset arithmetic. +// the 1D y-pass then uses gc[y*nx] with no extra offset arithmetic. inline void gradY(const float *v, double *vGrad, - double dy, + const double *y_, size_t nx, size_t ny, size_t nz) { - const double inv2dy = 1.0 / (2.0 * dy); const size_t slab = nx * ny; - for (int z = 0; z < (int)nz; ++z) { - for (int x = 0; x < (int)nx; ++x) { - const float *vc = v + (size_t)z * slab + x; - double *gc = vGrad + (size_t)z * slab + x; - gc[0] = (-3.0 * vc[0] + 4.0 * vc[nx] - vc[2*nx]) * inv2dy; - for (size_t y = 1; y < ny - 1; ++y) - gc[y * nx] = ((double)vc[(y+1)*nx] - (double)vc[(y-1)*nx]) * inv2dy; - gc[(ny-1) * nx] = ( 3.0 * vc[(ny-1)*nx] - 4.0 * vc[(ny-2)*nx] + vc[(ny-3)*nx]) * inv2dy; - } - } + for (int z = 0; z < (int)nz; ++z) + for (int x = 0; x < (int)nx; ++x) + grad1D_nu(v + (size_t)z*slab + x, vGrad + (size_t)z*slab + x, ny, nx, y_); } // Non-uniform z spacing: 1st-order one-sided at boundaries, 2nd-order central @@ -199,6 +228,63 @@ inline void gradZ(const float *w, } } +// --------------------------------------------------------------------------- +// vort_from_jacobians — compute vortical quantities from pre-computed Jacobian +// components. +// +// Inputs: u, v, w — float velocity components (any layout, length len) +// dux..dwz — Jacobian entries: d(u,v,w)/d(x,y,z), double, length len +// naming: d, e.g. duy = ∂u/∂y +// Outputs: vorticity, helicity, lambda2, qCriterion — float*, written in-place. +// Any output pointer may be null; null outputs are simply skipped. +// --------------------------------------------------------------------------- +inline void vort_from_jacobians(const float *u, + const float *v, + const float *w, + const double *dux, + const double *dvx, + const double *dwx, + const double *duy, + const double *dvy, + const double *dwy, + const double *duz, + const double *dvz, + const double *dwz, + float *vorticity, + float *helicity, + float *lambda2, + float *qCriterion, + size_t len) +{ + for (size_t i = 0; i < len; ++i) { + if (vorticity || helicity) { + const double omx = dwy[i] - dvz[i]; + const double omy = duz[i] - dwx[i]; + const double omz = dvx[i] - duy[i]; + const double omag = std::sqrt(omx * omx + omy * omy + omz * omz); + if (vorticity) + vorticity[i] = (float)omag; + if (helicity) { + const double ui = u[i], vi = v[i], wi = w[i]; + const double h = std::abs(omx * ui + omy * vi + omz * wi); + const double vmag = std::sqrt(ui * ui + vi * vi + wi * wi); + helicity[i] = (vmag > 0.0 && omag > 0.0) + ? (float)(h / (2.0 * vmag * omag)) + : 0.0f; + } + } + if (lambda2 || qCriterion) { + double J[3][3] = {{dux[i], duy[i], duz[i]}, + {dvx[i], dvy[i], dvz[i]}, + {dwx[i], dwy[i], dwz[i]}}; + if (lambda2) + lambda2[i] = (float)(-std::min(l2(J), 0.0)); + if (qCriterion) + qCriterion[i] = (float)std::max(q_crit(J), 0.0); + } + } +} + // --------------------------------------------------------------------------- // vort — compute vortical quantities from float velocity fields. // @@ -206,7 +292,6 @@ inline void gradZ(const float *w, // x_, y_, z_ — double coordinate arrays (length nx, ny, nz) // Outputs: vorticity, helicity, lambda2, qCriterion — float*, written in-place. // Any output pointer may be null; null outputs are simply skipped. -// If both lambda2 and qCriterion are null the eigenvalue loop is skipped. // // Internal gradient arrays are double-precision scratch space (9 × N doubles). // --------------------------------------------------------------------------- @@ -224,6 +309,9 @@ inline void vort(const float *u, size_t ny, size_t nz) { + if (!vorticity && !helicity && !lambda2 && !qCriterion) + return; + const size_t len = nx * ny * nz; // 9 double-precision gradient scratch arrays (shared across u/v/w) @@ -231,54 +319,35 @@ inline void vort(const float *u, std::vector duy(len), dvy(len), dwy(len); std::vector duz(len), dvz(len), dwz(len); - const double dx = x_[1] - x_[0]; - const double dy = y_[1] - y_[0]; - - gradX(u, dux.data(), dx, nx, ny, nz); - gradX(v, dvx.data(), dx, nx, ny, nz); - gradX(w, dwx.data(), dx, nx, ny, nz); + gradX(u, dux.data(), x_, nx, ny, nz); + gradX(v, dvx.data(), x_, nx, ny, nz); + gradX(w, dwx.data(), x_, nx, ny, nz); - gradY(u, duy.data(), dy, nx, ny, nz); - gradY(v, dvy.data(), dy, nx, ny, nz); - gradY(w, dwy.data(), dy, nx, ny, nz); + gradY(u, duy.data(), y_, nx, ny, nz); + gradY(v, dvy.data(), y_, nx, ny, nz); + gradY(w, dwy.data(), y_, nx, ny, nz); gradZ(u, duz.data(), z_, nx, ny, nz); gradZ(v, dvz.data(), z_, nx, ny, nz); gradZ(w, dwz.data(), z_, nx, ny, nz); - if (vorticity || helicity) { - for (size_t i = 0; i < len; ++i) { - const double omx = dwy[i] - dvz[i]; - const double omy = duz[i] - dwx[i]; - const double omz = dvx[i] - duy[i]; - const double omag = std::sqrt(omx * omx + omy * omy + omz * omz); - - if (vorticity) - vorticity[i] = (float)omag; - - if (helicity) { - const double ui = u[i], vi = v[i], wi = w[i]; - const double h = std::abs(omx * ui + omy * vi + omz * wi); - const double vmag = std::sqrt(ui * ui + vi * vi + wi * wi); - helicity[i] = (vmag > 0.0 && omag > 0.0) - ? (float)(h / (2.0 * vmag * omag)) - : 0.0f; - } - } - } - - // Lambda2 and Q-criterion (expensive: one eigensolve per voxel) - if (lambda2 || qCriterion) { - for (size_t i = 0; i < len; ++i) { - double J[3][3] = {{dux[i], duy[i], duz[i]}, - {dvx[i], dvy[i], dvz[i]}, - {dwx[i], dwy[i], dwz[i]}}; - if (lambda2) - lambda2[i] = (float)(-std::min(l2(J), 0.0)); - if (qCriterion) - qCriterion[i] = (float)std::max(q_crit(J), 0.0); - } - } + vort_from_jacobians(u, + v, + w, + dux.data(), + dvx.data(), + dwx.data(), + duy.data(), + dvy.data(), + dwy.data(), + duz.data(), + dvz.data(), + dwz.data(), + vorticity, + helicity, + lambda2, + qCriterion, + len); std::cout << "[vort] Vortical variables computed" << std::endl; } diff --git a/tsd/src/tsd/io/importers/import_VTU.cpp b/tsd/src/tsd/io/importers/import_VTU.cpp index 50e90e8fb..66052036b 100644 --- a/tsd/src/tsd/io/importers/import_VTU.cpp +++ b/tsd/src/tsd/io/importers/import_VTU.cpp @@ -257,7 +257,11 @@ static void createSurfaceFromGrid(Scene &scene, ("vtu_surface | " + std::string(filename)).c_str(), mesh, mat)); } -// Build an unstructured spatial field from volume cells only +// Build an unstructured spatial field from volume cells only. +// All point-data arrays are exposed as named SpatialFields sharing topology: +// 1-component → one field named after the array +// 3-component → three fields named {arr}_x, {arr}_y, {arr}_z +// Returns the first SpatialField created (used for the scene volume wrapper). static SpatialFieldRef createFieldFromVolumeCells( Scene &scene, vtkUnstructuredGrid *grid, const char *filepath) { @@ -276,11 +280,7 @@ static SpatialFieldRef createFieldFromVolumeCells( vtkIdType numPoints = grid->GetNumberOfPoints(); vtkIdType numVolumeCells = static_cast(volumeCellIndices.size()); - auto field = scene.createObject( - tokens::spatial_field::unstructured); - field->setName(fileOf(filepath).c_str()); - - // Vertices + // Build shared topology arrays auto vertexArray = scene.createArray(ANARI_FLOAT32_VEC3, numPoints); auto *vertexData = vertexArray->mapAs(); for (vtkIdType i = 0; i < numPoints; ++i) { @@ -288,9 +288,7 @@ static SpatialFieldRef createFieldFromVolumeCells( vertexData[i] = tsd::math::float3(pt[0], pt[1], pt[2]); } vertexArray->unmap(); - field->setParameterObject("vertex.position", *vertexArray); - // Cells std::vector connectivity; std::vector cellIndex; std::vector cellTypes; @@ -324,27 +322,84 @@ static SpatialFieldRef createFieldFromVolumeCells( auto indexArray = scene.createArray(ANARI_UINT32, connectivity.size()); indexArray->setData(connectivity.data()); - field->setParameterObject("index", *indexArray); auto cellIndexArray = scene.createArray(ANARI_UINT32, cellIndex.size()); cellIndexArray->setData(cellIndex.data()); - field->setParameterObject("cell.index", *cellIndexArray); auto cellTypesArray = scene.createArray(ANARI_UINT8, cellTypes.size()); cellTypesArray->setData(cellTypes.data()); - field->setParameterObject("cell.type", *cellTypesArray); - // Vertex data + // Helper: create a new unstructured SpatialField with shared topology + auto makeTopoField = [&](const std::string &name) -> SpatialFieldRef { + auto f = scene.createObject( + tokens::spatial_field::unstructured); + f->setName(name.c_str()); + f->setParameterObject("vertex.position", *vertexArray); + f->setParameterObject("index", *indexArray); + f->setParameterObject("cell.index", *cellIndexArray); + f->setParameterObject("cell.type", *cellTypesArray); + return f; + }; + + // Vertex data: expose all point arrays vtkPointData *pointData = grid->GetPointData(); - for (int i = 0; i < std::min(1, pointData->GetNumberOfArrays()); ++i) { + SpatialFieldRef firstField; + std::string baseName = fileOf(filepath); + + for (int i = 0; i < pointData->GetNumberOfArrays(); ++i) { vtkDataArray *array = pointData->GetArray(i); if (!array) continue; - auto a = makeFloatArray1D(scene, array, numPoints); - field->setParameterObject("vertex.data", *a); + int nComp = array->GetNumberOfComponents(); + std::string arrName = + (array->GetName() && array->GetName()[0] != '\0') ? array->GetName() + : baseName; + + if (nComp == 1) { + auto dataArr = scene.createArray(ANARI_FLOAT32, numPoints); + auto *buf = dataArr->mapAs(); + for (vtkIdType j = 0; j < numPoints; ++j) + buf[j] = static_cast(array->GetComponent(j, 0)); + dataArr->unmap(); + auto f = makeTopoField(arrName); + f->setParameterObject("vertex.data", *dataArr); + if (!firstField) + firstField = f; + } else if (nComp == 3) { + for (int c = 0; c < 3; ++c) { + const char *suffix = (c == 0) ? "_x" : (c == 1) ? "_y" : "_z"; + std::string compName = arrName + suffix; + auto dataArr = scene.createArray(ANARI_FLOAT32, numPoints); + auto *buf = dataArr->mapAs(); + for (vtkIdType j = 0; j < numPoints; ++j) + buf[j] = static_cast(array->GetComponent(j, c)); + dataArr->unmap(); + auto f = makeTopoField(compName); + f->setParameterObject("vertex.data", *dataArr); + if (!firstField) + firstField = f; + } + logStatus( + "[import_VTU] split 3-component array '%s' into '%s_x', '%s_y', " + "'%s_z'", + arrName.c_str(), + arrName.c_str(), + arrName.c_str(), + arrName.c_str()); + } else { + logWarning( + "[import_VTU] array '%s' has %d components (only 1 or 3 are " + "supported) -- skipping", + arrName.c_str(), + nComp); + } } - // Cell data + // Fallback: no point arrays → create a topology-only field + if (!firstField) + firstField = makeTopoField(baseName); + + // Cell data: attach first scalar array to firstField for rendering vtkCellData *cellData = grid->GetCellData(); for (int i = 0; i < std::min(1, cellData->GetNumberOfArrays()); ++i) { vtkDataArray *array = cellData->GetArray(i); @@ -361,10 +416,10 @@ static SpatialFieldRef createFieldFromVolumeCells( a->unmap(); filtered->unmap(); - field->setParameterObject("cell.data", *filtered); + firstField->setParameterObject("cell.data", *filtered); } - return field; + return firstField; } // Full-scene importer: surfaces + volumes diff --git a/tsd/src/tsd/io/procedural/computeVorticity.cpp b/tsd/src/tsd/io/procedural/computeVorticity.cpp index c14f138c9..34b00212c 100644 --- a/tsd/src/tsd/io/procedural/computeVorticity.cpp +++ b/tsd/src/tsd/io/procedural/computeVorticity.cpp @@ -13,6 +13,16 @@ // std #include #include +#if TSD_USE_VTK +#include +#include +#include +#include +#include +#include +#include +#include +#endif namespace tsd::io { @@ -172,6 +182,159 @@ static bool extractNanoVDB(Scene &scene, SpatialField *field, FieldData &out) return true; } +static bool extractStructuredRectilinear( + Scene &scene, SpatialField *field, FieldData &out) +{ + auto *p = field->parameter("data"); + if (!p || !anari::isArray(p->value().type())) { + logError( + "[computeVorticity] structuredRectilinear field has no 'data' array"); + return false; + } + + auto arr = scene.getObject(p->value().getAsObjectIndex()); + if (!arr) { + logError( + "[computeVorticity] structuredRectilinear 'data' array is invalid"); + return false; + } + + out.nx = arr->dim(0); + out.ny = arr->dim(1); + out.nz = arr->dim(2); + + if (out.nx < 2 || out.ny < 2 || out.nz < 2) { + logError("[computeVorticity] field dimensions must be >= 2 in each axis"); + return false; + } + + if (arr->elementType() != ANARI_FLOAT32) { + logError( + "[computeVorticity] structuredRectilinear field element type must be " + "ANARI_FLOAT32 (got %d)", + (int)arr->elementType()); + return false; + } + + out.ptr = arr->dataAs(); + + auto readCoords = [&](const char *name, + std::vector &dst, + size_t n) -> bool { + auto *cp = field->parameter(name); + if (!cp || !anari::isArray(cp->value().type())) { + logError( + "[computeVorticity] structuredRectilinear field missing '%s'", name); + return false; + } + auto ca = scene.getObject(cp->value().getAsObjectIndex()); + if (!ca) { + logError( + "[computeVorticity] structuredRectilinear '%s' array is invalid", + name); + return false; + } + dst.resize(n); + const float *src = ca->dataAs(); + for (size_t i = 0; i < n; ++i) + dst[i] = src[i]; + return true; + }; + + return readCoords("coordsX", out.x, out.nx) + && readCoords("coordsY", out.y, out.ny) + && readCoords("coordsZ", out.z, out.nz); +} + +static bool extractNanoVDBRectilinear( + Scene &scene, SpatialField *field, FieldData &out) +{ + auto *p = field->parameter("data"); + if (!p || !anari::isArray(p->value().type())) { + logError( + "[computeVorticity] nanovdbRectilinear field has no 'data' array"); + return false; + } + + auto arr = scene.getObject(p->value().getAsObjectIndex()); + if (!arr || !arr->data()) { + logError( + "[computeVorticity] nanovdbRectilinear 'data' array is invalid"); + return false; + } + + const uint8_t *rawData = static_cast(arr->data()); + const auto *meta = + reinterpret_cast(rawData); + + if (meta->gridType() != nanovdb::GridType::Float) { + logError( + "[computeVorticity] nanovdbRectilinear field must be float32 " + "(GridType::Float)"); + return false; + } + + const auto *grid = + reinterpret_cast *>(rawData); + auto bbox = grid->indexBBox(); + auto lo = bbox.min(); + auto hi = bbox.max(); + + out.nx = (size_t)(hi[0] - lo[0] + 1); + out.ny = (size_t)(hi[1] - lo[1] + 1); + out.nz = (size_t)(hi[2] - lo[2] + 1); + + if (out.nx < 2 || out.ny < 2 || out.nz < 2) { + logError( + "[computeVorticity] nanovdbRectilinear field dimensions must be >= 2"); + return false; + } + + auto readCoords = [&](const char *name, + std::vector &dst, + size_t n) -> bool { + auto *cp = field->parameter(name); + if (!cp || !anari::isArray(cp->value().type())) { + logError( + "[computeVorticity] nanovdbRectilinear field missing '%s'", name); + return false; + } + auto ca = scene.getObject(cp->value().getAsObjectIndex()); + if (!ca) { + logError( + "[computeVorticity] nanovdbRectilinear '%s' array is invalid", name); + return false; + } + dst.resize(n); + const float *src = ca->dataAs(); + for (size_t i = 0; i < n; ++i) + dst[i] = src[i]; + return true; + }; + + if (!readCoords("coordsX", out.x, out.nx) + || !readCoords("coordsY", out.y, out.ny) + || !readCoords("coordsZ", out.z, out.nz)) + return false; + + // Rasterize sparse grid into a dense float buffer + out.ownedData.resize(out.nx * out.ny * out.nz, 0.0f); + out.ptr = out.ownedData.data(); + auto acc = grid->getAccessor(); + + for (int k = lo[2]; k <= hi[2]; ++k) { + for (int j = lo[1]; j <= hi[1]; ++j) { + float *row = out.ownedData.data() + + (size_t)(k - lo[2]) * out.ny * out.nx + + (size_t)(j - lo[1]) * out.nx; + for (int i = lo[0]; i <= hi[0]; ++i) + row[i - lo[0]] = acc.getValue(nanovdb::Coord(i, j, k)); + } + } + + return true; +} + static bool extractFieldData(Scene &scene, SpatialField *field, FieldData &out) { if (!field) { @@ -181,12 +344,18 @@ static bool extractFieldData(Scene &scene, SpatialField *field, FieldData &out) if (field->subtype() == tokens::spatial_field::structuredRegular) { return extractStructuredRegular(scene, field, out); + } else if (field->subtype() == tokens::spatial_field::structuredRectilinear) { + return extractStructuredRectilinear(scene, field, out); } else if (field->subtype() == tokens::spatial_field::nanovdb) { return extractNanoVDB(scene, field, out); + } else if (field->subtype() == tokens::spatial_field::nanovdbRectilinear) { + return extractNanoVDBRectilinear(scene, field, out); } else { logError( "[computeVorticity] unsupported SpatialField subtype '%s'; " - "only structuredRegular and nanovdb are supported", + "only structuredRegular, structuredRectilinear, nanovdb, and " + "nanovdbRectilinear are supported (unstructured and amr require " + "mesh interpolation first)", field->subtype().c_str()); return false; } @@ -232,6 +401,265 @@ static VolumeRef wrapAsVolume(Scene &scene, return vol; } +// --------------------------------------------------------------------------- +// Unstructured path — only compiled when VTK is available +// --------------------------------------------------------------------------- + +#if TSD_USE_VTK + +// Create an unstructured SpatialField + transferFunction1D Volume by copying +// topology parameters from meshSource and supplying a new vertex.data array. +static VolumeRef wrapAsUnstructuredVolume(Scene &scene, + const std::string &name, + SpatialField *meshSource, + ArrayRef &dataArr, + LayerNodeRef location) +{ + auto field = scene.createObject( + tokens::spatial_field::unstructured); + field->setName(name.c_str()); + + // Share topology parameter objects from meshSource + for (const char *pname : + {"vertex.position", "index", "cell.index", "cell.type"}) { + auto *p = meshSource->parameter(pname); + if (!p || !anari::isArray(p->value().type())) + continue; + auto arr = scene.getObject(p->value().getAsObjectIndex()); + if (arr) + field->setParameterObject(pname, *arr); + } + + field->setParameterObject("vertex.data", *dataArr); + + float2 valueRange = field->computeValueRange(); + + auto tx = scene.insertChildTransformNode( + location ? location : scene.defaultLayer()->root()); + + auto [inst, vol] = + scene.insertNewChildObjectNode(tx, tokens::volume::transferFunction1D); + vol->setName(name.c_str()); + vol->setParameterObject("value", *field); + vol->setParameter("valueRange", ANARI_FLOAT32_BOX1, &valueRange); + + auto colorArr = scene.createArray(ANARI_FLOAT32_VEC4, 256); + colorArr->setData(makeDefaultColorMap(256).data()); + vol->setParameterObject("color", *colorArr); + + return vol; +} + +static VorticityResult computeVorticityUnstructured(Scene &scene, + SpatialField *u, + SpatialField *v, + SpatialField *w, + LayerNodeRef location, + VorticityOptions opts) +{ + VorticityResult result; + + // Step 1: numPoints from u's vertex.position array + auto *posParam = u->parameter("vertex.position"); + if (!posParam || !anari::isArray(posParam->value().type())) { + logError( + "[computeVorticity] unstructured field missing 'vertex.position'"); + return result; + } + auto posArr = scene.getObject(posParam->value().getAsObjectIndex()); + if (!posArr) { + logError( + "[computeVorticity] unstructured field 'vertex.position' is invalid"); + return result; + } + const size_t numPoints = posArr->dim(0); + + // Validate v and w have matching point count + auto getNumPoints = [&](SpatialField *field) -> size_t { + auto *p = field->parameter("vertex.position"); + if (!p || !anari::isArray(p->value().type())) + return 0; + auto a = scene.getObject(p->value().getAsObjectIndex()); + return a ? a->dim(0) : 0; + }; + if (getNumPoints(v) != numPoints || getNumPoints(w) != numPoints) { + logError( + "[computeVorticity] U/V/W unstructured fields have mismatched point " + "counts"); + return result; + } + + // Step 2: scalar velocity data from vertex.data + auto getDataPtr = [&](SpatialField *field) -> const float * { + auto *p = field->parameter("vertex.data"); + if (!p || !anari::isArray(p->value().type())) + return nullptr; + auto a = scene.getObject(p->value().getAsObjectIndex()); + return a ? a->dataAs() : nullptr; + }; + const float *uPtr = getDataPtr(u); + const float *vPtr = getDataPtr(v); + const float *wPtr = getDataPtr(w); + if (!uPtr || !vPtr || !wPtr) { + logError( + "[computeVorticity] unstructured fields missing 'vertex.data' scalar"); + return result; + } + + // Step 3: get topology arrays + auto *indexParam = u->parameter("index"); + auto *cellIdxParam = u->parameter("cell.index"); + auto *cellTypeParam = u->parameter("cell.type"); + if (!indexParam || !cellIdxParam || !cellTypeParam + || !anari::isArray(indexParam->value().type()) + || !anari::isArray(cellIdxParam->value().type()) + || !anari::isArray(cellTypeParam->value().type())) { + logError( + "[computeVorticity] unstructured field missing topology parameters"); + return result; + } + auto connArr = scene.getObject(indexParam->value().getAsObjectIndex()); + auto cellIdxArr = + scene.getObject(cellIdxParam->value().getAsObjectIndex()); + auto cellTypeArr = + scene.getObject(cellTypeParam->value().getAsObjectIndex()); + + const uint32_t *connectivity = connArr->dataAs(); + const uint32_t *cellIndex = cellIdxArr->dataAs(); + const uint8_t *cellTypeData = cellTypeArr->dataAs(); + const size_t numCells = cellIdxArr->dim(0); + const size_t connSize = connArr->dim(0); + + logStatus( + "[computeVorticity] computing vortical quantities on %zu-point " + "unstructured grid...", + numPoints); + + // Step 4: reconstruct vtkUnstructuredGrid + auto vgrid = vtkSmartPointer::New(); + + auto points = vtkSmartPointer::New(); + points->SetDataTypeToFloat(); + points->SetNumberOfPoints(static_cast(numPoints)); + const float *rawPos = posArr->dataAs(); + for (size_t i = 0; i < numPoints; ++i) + points->SetPoint(static_cast(i), + rawPos[i * 3], + rawPos[i * 3 + 1], + rawPos[i * 3 + 2]); + vgrid->SetPoints(points); + + vgrid->Allocate(static_cast(numCells)); + for (size_t ci = 0; ci < numCells; ++ci) { + size_t start = cellIndex[ci]; + size_t end = (ci + 1 < numCells) ? cellIndex[ci + 1] : connSize; + size_t npts = end - start; + std::vector ptIds(npts); + for (size_t j = 0; j < npts; ++j) + ptIds[j] = static_cast(connectivity[start + j]); + vgrid->InsertNextCell( + static_cast(cellTypeData[ci]), static_cast(npts), ptIds.data()); + } + + // Add interleaved velocity array (u,v,w per point) + auto velArr = vtkSmartPointer::New(); + velArr->SetName("velocity"); + velArr->SetNumberOfComponents(3); + velArr->SetNumberOfTuples(static_cast(numPoints)); + for (size_t i = 0; i < numPoints; ++i) { + float tuple[3] = {uPtr[i], vPtr[i], wPtr[i]}; + velArr->SetTuple(static_cast(i), tuple); + } + vgrid->GetPointData()->AddArray(velArr); + + // Step 5: run vtkGradientFilter to get the 9-component Jacobian + auto gradFilter = vtkSmartPointer::New(); + gradFilter->SetInputData(vgrid); + gradFilter->SetInputArrayToProcess( + 0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, "velocity"); + gradFilter->SetResultArrayName("gradients"); + gradFilter->SetComputeGradient(1); + gradFilter->Update(); + + auto *gradOutput = + vtkUnstructuredGrid::SafeDownCast(gradFilter->GetOutput()); + if (!gradOutput) { + logError("[computeVorticity] vtkGradientFilter failed"); + return result; + } + vtkDataArray *gradData = + gradOutput->GetPointData()->GetArray("gradients"); + if (!gradData || gradData->GetNumberOfComponents() != 9) { + logError( + "[computeVorticity] gradient filter output has unexpected format"); + return result; + } + + // Step 6: unpack Jacobian into 9 double vectors + // VTK ordering: d(comp0)/dx, d(comp0)/dy, d(comp0)/dz, + // d(comp1)/dx, d(comp1)/dy, d(comp1)/dz, + // d(comp2)/dx, d(comp2)/dy, d(comp2)/dz + // → dux, duy, duz, dvx, dvy, dvz, dwx, dwy, dwz + std::vector dux(numPoints), duy(numPoints), duz(numPoints); + std::vector dvx(numPoints), dvy(numPoints), dvz(numPoints); + std::vector dwx(numPoints), dwy(numPoints), dwz(numPoints); + for (size_t i = 0; i < numPoints; ++i) { + double t[9]; + gradData->GetTuple(static_cast(i), t); + dux[i] = t[0]; duy[i] = t[1]; duz[i] = t[2]; + dvx[i] = t[3]; dvy[i] = t[4]; dvz[i] = t[5]; + dwx[i] = t[6]; dwy[i] = t[7]; dwz[i] = t[8]; + } + + // Step 7: allocate output arrays and compute + auto makeOutBuf = [&](bool enabled) -> std::pair { + if (!enabled) + return {}; + auto arr = scene.createArray(ANARI_FLOAT32, numPoints); + return {arr, arr->mapAs()}; + }; + auto [lambda2Arr, lambda2Out] = makeOutBuf(opts.lambda2); + auto [qCritArr, qCritOut] = makeOutBuf(opts.qCriterion); + auto [vorticityArr, vorticityOut] = makeOutBuf(opts.vorticity); + auto [helicityArr, helicityOut] = makeOutBuf(opts.helicity); + + vort_from_jacobians(uPtr, + vPtr, + wPtr, + dux.data(), dvx.data(), dwx.data(), + duy.data(), dvy.data(), dwy.data(), + duz.data(), dvz.data(), dwz.data(), + vorticityOut, helicityOut, lambda2Out, qCritOut, + numPoints); + + // Step 8: unmap + if (lambda2Arr) lambda2Arr->unmap(); + if (qCritArr) qCritArr->unmap(); + if (vorticityArr) vorticityArr->unmap(); + if (helicityArr) helicityArr->unmap(); + + logStatus("[computeVorticity] creating output volumes..."); + + // Step 9: wrap each result as an unstructured volume + if (lambda2Arr) + result.lambda2 = + wrapAsUnstructuredVolume(scene, "lambda2", u, lambda2Arr, location); + if (qCritArr) + result.qCriterion = + wrapAsUnstructuredVolume(scene, "q_criterion", u, qCritArr, location); + if (vorticityArr) + result.vorticity = + wrapAsUnstructuredVolume(scene, "vorticity", u, vorticityArr, location); + if (helicityArr) + result.helicity = + wrapAsUnstructuredVolume(scene, "helicity", u, helicityArr, location); + + logStatus("[computeVorticity] done."); + return result; +} + +#endif // TSD_USE_VTK + // --------------------------------------------------------------------------- // Public API // --------------------------------------------------------------------------- @@ -245,6 +673,11 @@ VorticityResult computeVorticity(Scene &scene, { VorticityResult result; +#if TSD_USE_VTK + if (u && u->subtype() == tokens::spatial_field::unstructured) + return computeVorticityUnstructured(scene, u, v, w, location, opts); +#endif + logStatus("[computeVorticity] extracting field data..."); FieldData uData, vData, wData; @@ -275,26 +708,16 @@ VorticityResult computeVorticity(Scene &scene, // Pre-allocate one ANARI float32 array per selected output and map it so // vort() can write directly into the final storage — no intermediate buffers. - ArrayRef lambda2Arr, qCritArr, vorticityArr, helicityArr; - float *lambda2Out = nullptr, *qCritOut = nullptr; - float *vorticityOut = nullptr, *helicityOut = nullptr; - - if (opts.lambda2) { - lambda2Arr = scene.createArray(ANARI_FLOAT32, nx, ny, nz); - lambda2Out = lambda2Arr->mapAs(); - } - if (opts.qCriterion) { - qCritArr = scene.createArray(ANARI_FLOAT32, nx, ny, nz); - qCritOut = qCritArr->mapAs(); - } - if (opts.vorticity) { - vorticityArr = scene.createArray(ANARI_FLOAT32, nx, ny, nz); - vorticityOut = vorticityArr->mapAs(); - } - if (opts.helicity) { - helicityArr = scene.createArray(ANARI_FLOAT32, nx, ny, nz); - helicityOut = helicityArr->mapAs(); - } + auto makeOutBuf = [&](bool enabled) -> std::pair { + if (!enabled) + return {}; + auto arr = scene.createArray(ANARI_FLOAT32, nx, ny, nz); + return {arr, arr->mapAs()}; + }; + auto [lambda2Arr, lambda2Out] = makeOutBuf(opts.lambda2); + auto [qCritArr, qCritOut] = makeOutBuf(opts.qCriterion); + auto [vorticityArr, vorticityOut] = makeOutBuf(opts.vorticity); + auto [helicityArr, helicityOut] = makeOutBuf(opts.helicity); // Compute — vort() reads float* velocity directly, writes float* outputs. // Null output pointers are silently skipped inside vort(). From 6d3c62491837eb4ce29ce26d3e5d009a4cbcfb80 Mon Sep 17 00:00:00 2001 From: Andrea Paris Date: Fri, 27 Feb 2026 16:24:06 +0100 Subject: [PATCH 05/12] VTU volume TF fix --- tsd/src/tsd/io/importers/import_VTU.cpp | 3 + .../tsd/io/procedural/computeVorticity.cpp | 156 +++++++++++------- 2 files changed, 98 insertions(+), 61 deletions(-) diff --git a/tsd/src/tsd/io/importers/import_VTU.cpp b/tsd/src/tsd/io/importers/import_VTU.cpp index 66052036b..40de88214 100644 --- a/tsd/src/tsd/io/importers/import_VTU.cpp +++ b/tsd/src/tsd/io/importers/import_VTU.cpp @@ -466,6 +466,9 @@ void import_VTU(Scene &scene, const char *filepath, LayerNodeRef location) volume->setName(fileOf(filepath).c_str()); volume->setParameterObject("value", *field); volume->setParameter("valueRange", ANARI_FLOAT32_BOX1, &valueRange); + auto colorArr = scene.createArray(ANARI_FLOAT32_VEC4, 256); + colorArr->setData(tsd::core::makeDefaultColorMap(256).data()); + volume->setParameterObject("color", *colorArr); } } } diff --git a/tsd/src/tsd/io/procedural/computeVorticity.cpp b/tsd/src/tsd/io/procedural/computeVorticity.cpp index 34b00212c..09ad3d913 100644 --- a/tsd/src/tsd/io/procedural/computeVorticity.cpp +++ b/tsd/src/tsd/io/procedural/computeVorticity.cpp @@ -18,8 +18,10 @@ #include #include #include +#include #include #include +#include #include #include #endif @@ -407,48 +409,6 @@ static VolumeRef wrapAsVolume(Scene &scene, #if TSD_USE_VTK -// Create an unstructured SpatialField + transferFunction1D Volume by copying -// topology parameters from meshSource and supplying a new vertex.data array. -static VolumeRef wrapAsUnstructuredVolume(Scene &scene, - const std::string &name, - SpatialField *meshSource, - ArrayRef &dataArr, - LayerNodeRef location) -{ - auto field = scene.createObject( - tokens::spatial_field::unstructured); - field->setName(name.c_str()); - - // Share topology parameter objects from meshSource - for (const char *pname : - {"vertex.position", "index", "cell.index", "cell.type"}) { - auto *p = meshSource->parameter(pname); - if (!p || !anari::isArray(p->value().type())) - continue; - auto arr = scene.getObject(p->value().getAsObjectIndex()); - if (arr) - field->setParameterObject(pname, *arr); - } - - field->setParameterObject("vertex.data", *dataArr); - - float2 valueRange = field->computeValueRange(); - - auto tx = scene.insertChildTransformNode( - location ? location : scene.defaultLayer()->root()); - - auto [inst, vol] = - scene.insertNewChildObjectNode(tx, tokens::volume::transferFunction1D); - vol->setName(name.c_str()); - vol->setParameterObject("value", *field); - vol->setParameter("valueRange", ANARI_FLOAT32_BOX1, &valueRange); - - auto colorArr = scene.createArray(ANARI_FLOAT32_VEC4, 256); - colorArr->setData(makeDefaultColorMap(256).data()); - vol->setParameterObject("color", *colorArr); - - return vol; -} static VorticityResult computeVorticityUnstructured(Scene &scene, SpatialField *u, @@ -541,12 +501,10 @@ static VorticityResult computeVorticityUnstructured(Scene &scene, auto points = vtkSmartPointer::New(); points->SetDataTypeToFloat(); points->SetNumberOfPoints(static_cast(numPoints)); - const float *rawPos = posArr->dataAs(); + const auto *rawPos = posArr->dataAs(); for (size_t i = 0; i < numPoints; ++i) - points->SetPoint(static_cast(i), - rawPos[i * 3], - rawPos[i * 3 + 1], - rawPos[i * 3 + 2]); + points->SetPoint( + static_cast(i), rawPos[i].x, rawPos[i].y, rawPos[i].z); vgrid->SetPoints(points); vgrid->Allocate(static_cast(numCells)); @@ -638,21 +596,97 @@ static VorticityResult computeVorticityUnstructured(Scene &scene, if (vorticityArr) vorticityArr->unmap(); if (helicityArr) helicityArr->unmap(); - logStatus("[computeVorticity] creating output volumes..."); + // Step 9: Resample onto a structured regular grid and create volumes. + // VisRTX does not support the "unstructured" SpatialField subtype, so we + // interpolate the per-point output values onto an axis-aligned regular grid + // using vtkProbeFilter, then wrap as structuredRegular which VisRTX does + // support. + + // 9a. Add vortical output arrays back to vgrid as VTK point data. + auto addToGrid = [&](const char *name, ArrayRef &arr, bool enabled) { + if (!enabled || !arr) + return; + const float *src = arr->dataAs(); // ANARI_FLOAT32, safe + auto vtkarr = vtkSmartPointer::New(); + vtkarr->SetName(name); + vtkarr->SetNumberOfComponents(1); + vtkarr->SetNumberOfTuples(static_cast(numPoints)); + for (size_t i = 0; i < numPoints; ++i) + vtkarr->SetValue(static_cast(i), src[i]); + vgrid->GetPointData()->AddArray(vtkarr); + }; + addToGrid("lambda2", lambda2Arr, opts.lambda2); + addToGrid("q_criterion", qCritArr, opts.qCriterion); + addToGrid("vorticity", vorticityArr, opts.vorticity); + addToGrid("helicity", helicityArr, opts.helicity); + + // 9b. Compute output grid resolution proportional to bounding box extents. + double bounds[6]; + vgrid->GetBounds(bounds); + const double bx = bounds[1] - bounds[0]; + const double by = bounds[3] - bounds[2]; + const double bz = bounds[5] - bounds[4]; + const double maxB = std::max({bx, by, bz}); + if (maxB < 1e-10) { + logError("[computeVorticity] degenerate mesh bounds, aborting resample"); + return result; + } + const int MAX_RES = 64; + const size_t resX = std::max(size_t(2), size_t(std::round(MAX_RES * bx / maxB))); + const size_t resY = std::max(size_t(2), size_t(std::round(MAX_RES * by / maxB))); + const size_t resZ = std::max(size_t(2), size_t(std::round(MAX_RES * bz / maxB))); + + logStatus("[computeVorticity] resampling to %zu×%zu×%zu structured grid...", + resX, resY, resZ); + + auto imData = vtkSmartPointer::New(); + imData->SetDimensions(static_cast(resX), + static_cast(resY), + static_cast(resZ)); + imData->SetOrigin(bounds[0], bounds[2], bounds[4]); + imData->SetSpacing( + bx / (resX - 1), by / (resY - 1), bz / (resZ - 1)); + + auto probe = vtkSmartPointer::New(); + probe->SetSourceData(vgrid); + probe->SetInputData(imData); + probe->Update(); + + auto *probed = vtkImageData::SafeDownCast(probe->GetOutput()); + if (!probed) { + logError("[computeVorticity] vtkProbeFilter resampling failed"); + return result; + } - // Step 9: wrap each result as an unstructured volume - if (lambda2Arr) - result.lambda2 = - wrapAsUnstructuredVolume(scene, "lambda2", u, lambda2Arr, location); - if (qCritArr) - result.qCriterion = - wrapAsUnstructuredVolume(scene, "q_criterion", u, qCritArr, location); - if (vorticityArr) - result.vorticity = - wrapAsUnstructuredVolume(scene, "vorticity", u, vorticityArr, location); - if (helicityArr) - result.helicity = - wrapAsUnstructuredVolume(scene, "helicity", u, helicityArr, location); + // 9c. Extract resampled arrays and wrap as structuredRegular volumes. + const size_t resTotal = resX * resY * resZ; + const math::float3 origin{ + (float)bounds[0], (float)bounds[2], (float)bounds[4]}; + const math::float3 spacing{ + (float)(bx / (resX - 1)), + (float)(by / (resY - 1)), + (float)(bz / (resZ - 1))}; + + auto wrapResampled = [&](const char *name, bool enabled) -> VolumeRef { + if (!enabled) + return {}; + vtkDataArray *arr = probed->GetPointData()->GetArray(name); + if (!arr) + return {}; + auto outArr = scene.createArray(ANARI_FLOAT32, resX, resY, resZ); + auto *buf = outArr->mapAs(); + for (size_t i = 0; i < resTotal; ++i) + buf[i] = static_cast(arr->GetTuple1(static_cast(i))); + outArr->unmap(); + return wrapAsVolume( + scene, name, outArr, resX, resY, resZ, origin, spacing, location); + }; + + logStatus("[computeVorticity] creating output volumes..."); + result.lambda2 = wrapResampled("lambda2", opts.lambda2); + result.qCriterion = wrapResampled("q_criterion", opts.qCriterion); + result.vorticity = wrapResampled("vorticity", opts.vorticity); + result.helicity = wrapResampled("helicity", opts.helicity); logStatus("[computeVorticity] done."); return result; From 31a78526710304272e0ffaa98fec9ea6e6530404 Mon Sep 17 00:00:00 2001 From: Andrea Paris Date: Fri, 27 Feb 2026 17:25:55 +0100 Subject: [PATCH 06/12] lambda2 cleanup --- tsd/src/tsd/core/algorithms/vort.h | 129 +++++++---------------------- 1 file changed, 28 insertions(+), 101 deletions(-) diff --git a/tsd/src/tsd/core/algorithms/vort.h b/tsd/src/tsd/core/algorithms/vort.h index bcc9bebed..fe64a3bd8 100644 --- a/tsd/src/tsd/core/algorithms/vort.h +++ b/tsd/src/tsd/core/algorithms/vort.h @@ -8,66 +8,45 @@ #include #include -// 3x3 symmetric eigenvalue solver — analytical (Cardano/trigonometric method) -// Returns eigenvalues sorted largest→smallest in eigs[0..2] -static void eigen3sym(double A[3][3], double eigs[3]) -{ - double p1 = A[0][1] * A[0][1] + A[0][2] * A[0][2] + A[1][2] * A[1][2]; - if (p1 == 0.0) { - // Diagonal matrix — eigenvalues are the diagonal entries - eigs[0] = A[0][0]; - eigs[1] = A[1][1]; - eigs[2] = A[2][2]; - if (eigs[0] < eigs[1]) - std::swap(eigs[0], eigs[1]); - if (eigs[1] < eigs[2]) - std::swap(eigs[1], eigs[2]); - if (eigs[0] < eigs[1]) - std::swap(eigs[0], eigs[1]); - return; - } - double q = (A[0][0] + A[1][1] + A[2][2]) / 3.0; - double p2 = (A[0][0] - q) * (A[0][0] - q) + (A[1][1] - q) * (A[1][1] - q) - + (A[2][2] - q) * (A[2][2] - q) + 2 * p1; - double p = std::sqrt(p2 / 6.0); - double B[3][3]; - for (int i = 0; i < 3; i++) - for (int j = 0; j < 3; j++) - B[i][j] = (A[i][j] - (i == j ? q : 0.0)) / p; - double r = (B[0][0] * (B[1][1] * B[2][2] - B[1][2] * B[2][1]) - - B[0][1] * (B[1][0] * B[2][2] - B[1][2] * B[2][0]) - + B[0][2] * (B[1][0] * B[2][1] - B[1][1] * B[2][0])) - / 2.0; - double phi; - if (r <= -1.0) - phi = M_PI / 3.0; - else if (r >= 1.0) - phi = 0.0; - else - phi = std::acos(r) / 3.0; - eigs[0] = q + 2 * p * std::cos(phi); - eigs[2] = q + 2 * p * std::cos(phi + 2.0943951023931953); // phi + 2*pi/3 - eigs[1] = 3 * q - eigs[0] - eigs[2]; -} - -// lambda2: second (middle) eigenvalue of S^2 + O^2 -// J is the velocity gradient Jacobian; uses no external dependencies +// lambda2: middle eigenvalue of M = S^2 + O^2, where S and O are the +// symmetric and skew-symmetric parts of the velocity gradient Jacobian J. +// Uses the analytical trigonometric (Cardano) formula for 3x3 symmetric matrices: +// eig_k = q + 2p * cos(phi + 2*pi*k/3), k = 0,1,2 +// where q = tr(M)/3, p = ||M - q*I||_F / sqrt(6), phi = acos(det(M-q*I)/(2p^3)) / 3. +// The middle eigenvalue is recovered via the trace identity: eig1 = 3q - eig0 - eig2. static double l2(double J[3][3]) { double S[3][3], O[3][3]; for (int i = 0; i < 3; i++) for (int j = 0; j < 3; j++) { S[i][j] = 0.5 * (J[i][j] + J[j][i]); - O[i][j] = S[i][j] - J[j][i]; // O = S - J^T = 0.5*(J - J^T) + O[i][j] = 0.5 * (J[i][j] - J[j][i]); } double M[3][3] = {}; for (int i = 0; i < 3; i++) for (int j = 0; j < 3; j++) for (int k = 0; k < 3; k++) M[i][j] += S[i][k] * S[k][j] + O[i][k] * O[k][j]; - double eigs[3]; - eigen3sym(M, eigs); // sorted largest→smallest; eigs[1] is the middle - return eigs[1]; + + const double q = (M[0][0] + M[1][1] + M[2][2]) / 3.0; + const double a = M[0][0]-q, d = M[1][1]-q, f = M[2][2]-q; + const double p1 = M[0][1]*M[0][1] + M[0][2]*M[0][2] + M[1][2]*M[1][2]; + if (p1 == 0.0) { + // Diagonal — middle eigenvalue by inspection + double e[3] = {a, d, f}; + if (e[0] < e[1]) std::swap(e[0], e[1]); + if (e[1] < e[2]) std::swap(e[1], e[2]); + if (e[0] < e[1]) std::swap(e[0], e[1]); + return q + e[1]; + } + const double p = std::sqrt((a*a + d*d + f*f + 2.0*p1) / 6.0); + const double r = (a*(d*f - M[1][2]*M[1][2]) + - M[0][1]*(M[0][1]*f - M[1][2]*M[0][2]) + + M[0][2]*(M[0][1]*M[1][2] - d*M[0][2])) / (2.0*p*p*p); + const double phi = std::acos(std::max(-1.0, std::min(1.0, r))) / 3.0; + const double e0 = q + 2.0*p*std::cos(phi); + const double e2 = q + 2.0*p*std::cos(phi + 2.0943951023931953); // phi + 2*pi/3 + return 3.0*q - e0 - e2; // middle eigenvalue via trace identity } // Q-criterion: 0.5 * (||O||_F^2 - ||S||_F^2) @@ -77,7 +56,7 @@ static double q_crit(double J[3][3]) for (int i = 0; i < 3; i++) for (int j = 0; j < 3; j++) { S[i][j] = 0.5 * (J[i][j] + J[j][i]); - O[i][j] = S[i][j] - J[j][i]; + O[i][j] = 0.5 * (J[i][j] - J[j][i]); } double trO2 = 0, trS2 = 0; for (int i = 0; i < 3; i++) @@ -124,58 +103,6 @@ static void grad1D_nu( // Boundary accuracy: 1st-order one-sided. Interior: 2nd-order central. // --------------------------------------------------------------------------- -// Forward declarations -inline void gradX(const float *u, - double *uGrad, - const double *x_, - size_t nx, - size_t ny, - size_t nz); -inline void gradY(const float *v, - double *vGrad, - const double *y_, - size_t nx, - size_t ny, - size_t nz); -inline void gradZ(const float *w, - double *wGrad, - const double *z_, - size_t nx, - size_t ny, - size_t nz); -inline void vort_from_jacobians(const float *u, - const float *v, - const float *w, - const double *dux, - const double *dvx, - const double *dwx, - const double *duy, - const double *dvy, - const double *dwy, - const double *duz, - const double *dvz, - const double *dwz, - float *vorticity, - float *helicity, - float *lambda2, - float *qCriterion, - size_t len); -inline void vort(const float *u, - const float *v, - const float *w, - const double *x_, - const double *y_, - const double *z_, - float *vorticity, - float *helicity, - float *lambda2, - float *qCriterion, - size_t nx, - size_t ny, - size_t nz); - -// --------------------------------------------------------------------------- - // x is the fastest index, so every (z,y) pair is a contiguous row of nx // floats — collapse both outer loops into a single row index. inline void gradX(const float *u, From 0a7130c8fc34a09a79226cf9ca814ec086c6616d Mon Sep 17 00:00:00 2001 From: Andrea Paris Date: Fri, 27 Feb 2026 19:44:32 +0100 Subject: [PATCH 07/12] lambda2 cleanup --- tsd/src/tsd/core/algorithms/vort.h | 74 ++++++++++++------------------ 1 file changed, 30 insertions(+), 44 deletions(-) diff --git a/tsd/src/tsd/core/algorithms/vort.h b/tsd/src/tsd/core/algorithms/vort.h index fe64a3bd8..234432ad9 100644 --- a/tsd/src/tsd/core/algorithms/vort.h +++ b/tsd/src/tsd/core/algorithms/vort.h @@ -8,29 +8,25 @@ #include #include -// lambda2: middle eigenvalue of M = S^2 + O^2, where S and O are the -// symmetric and skew-symmetric parts of the velocity gradient Jacobian J. -// Uses the analytical trigonometric (Cardano) formula for 3x3 symmetric matrices: -// eig_k = q + 2p * cos(phi + 2*pi*k/3), k = 0,1,2 -// where q = tr(M)/3, p = ||M - q*I||_F / sqrt(6), phi = acos(det(M-q*I)/(2p^3)) / 3. -// The middle eigenvalue is recovered via the trace identity: eig1 = 3q - eig0 - eig2. -static double l2(double J[3][3]) +// lambda2: middle eigenvalue of M = S^2 + O^2, where S=(J+J^T)/2, O=(J-J^T)/2. +// M is symmetric with 6 unique entries computed as M = (J^2 + (J^T)^2) / 2. +// Eigenvalues via the trigonometric (Cardano) formula: eig_k = q + 2p*cos(phi + 2*pi*k/3). +// Middle eigenvalue recovered from the trace identity: eig1 = 3q - eig0 - eig2. +static double l2(double j00, double j01, double j02, + double j10, double j11, double j12, + double j20, double j21, double j22) { - double S[3][3], O[3][3]; - for (int i = 0; i < 3; i++) - for (int j = 0; j < 3; j++) { - S[i][j] = 0.5 * (J[i][j] + J[j][i]); - O[i][j] = 0.5 * (J[i][j] - J[j][i]); - } - double M[3][3] = {}; - for (int i = 0; i < 3; i++) - for (int j = 0; j < 3; j++) - for (int k = 0; k < 3; k++) - M[i][j] += S[i][k] * S[k][j] + O[i][k] * O[k][j]; + // 6 unique entries of symmetric M = (J^2 + (J^T)^2) / 2 + const double m00 = j00*j00 + j01*j10 + j02*j20; + const double m11 = j10*j01 + j11*j11 + j12*j21; + const double m22 = j20*j02 + j21*j12 + j22*j22; + const double m01 = 0.5*(j00*j01 + j01*j11 + j02*j21 + j00*j10 + j10*j11 + j20*j12); + const double m02 = 0.5*(j00*j02 + j01*j12 + j02*j22 + j00*j20 + j10*j21 + j20*j22); + const double m12 = 0.5*(j10*j02 + j11*j12 + j12*j22 + j01*j20 + j11*j21 + j21*j22); - const double q = (M[0][0] + M[1][1] + M[2][2]) / 3.0; - const double a = M[0][0]-q, d = M[1][1]-q, f = M[2][2]-q; - const double p1 = M[0][1]*M[0][1] + M[0][2]*M[0][2] + M[1][2]*M[1][2]; + const double q = (m00 + m11 + m22) / 3.0; + const double a = m00-q, d = m11-q, f = m22-q; + const double p1 = m01*m01 + m02*m02 + m12*m12; if (p1 == 0.0) { // Diagonal — middle eigenvalue by inspection double e[3] = {a, d, f}; @@ -40,31 +36,20 @@ static double l2(double J[3][3]) return q + e[1]; } const double p = std::sqrt((a*a + d*d + f*f + 2.0*p1) / 6.0); - const double r = (a*(d*f - M[1][2]*M[1][2]) - - M[0][1]*(M[0][1]*f - M[1][2]*M[0][2]) - + M[0][2]*(M[0][1]*M[1][2] - d*M[0][2])) / (2.0*p*p*p); + const double r = (a*(d*f - m12*m12) - m01*(m01*f - m12*m02) + m02*(m01*m12 - d*m02)) + / (2.0*p*p*p); const double phi = std::acos(std::max(-1.0, std::min(1.0, r))) / 3.0; const double e0 = q + 2.0*p*std::cos(phi); const double e2 = q + 2.0*p*std::cos(phi + 2.0943951023931953); // phi + 2*pi/3 return 3.0*q - e0 - e2; // middle eigenvalue via trace identity } -// Q-criterion: 0.5 * (||O||_F^2 - ||S||_F^2) -static double q_crit(double J[3][3]) +// Q-criterion: 0.5*(||O||^2 - ||S||^2) = -0.5*tr(J^2) = -0.5*(j00^2+j11^2+j22^2 + 2*(j01*j10+j02*j20+j12*j21)) +static double q_crit(double j00, double j01, double j02, + double j10, double j11, double j12, + double j20, double j21, double j22) { - double S[3][3], O[3][3]; - for (int i = 0; i < 3; i++) - for (int j = 0; j < 3; j++) { - S[i][j] = 0.5 * (J[i][j] + J[j][i]); - O[i][j] = 0.5 * (J[i][j] - J[j][i]); - } - double trO2 = 0, trS2 = 0; - for (int i = 0; i < 3; i++) - for (int j = 0; j < 3; j++) { - trO2 += O[i][j] * O[i][j]; - trS2 += S[i][j] * S[i][j]; - } - return 0.5 * (trO2 - trS2); + return -0.5*(j00*j00 + j11*j11 + j22*j22 + 2.0*(j01*j10 + j02*j20 + j12*j21)); } // --------------------------------------------------------------------------- @@ -201,13 +186,14 @@ inline void vort_from_jacobians(const float *u, } } if (lambda2 || qCriterion) { - double J[3][3] = {{dux[i], duy[i], duz[i]}, - {dvx[i], dvy[i], dvz[i]}, - {dwx[i], dwy[i], dwz[i]}}; if (lambda2) - lambda2[i] = (float)(-std::min(l2(J), 0.0)); + lambda2[i] = (float)(-std::min(l2(dux[i], duy[i], duz[i], + dvx[i], dvy[i], dvz[i], + dwx[i], dwy[i], dwz[i]), 0.0)); if (qCriterion) - qCriterion[i] = (float)std::max(q_crit(J), 0.0); + qCriterion[i] = (float)std::max(q_crit(dux[i], duy[i], duz[i], + dvx[i], dvy[i], dvz[i], + dwx[i], dwy[i], dwz[i]), 0.0); } } } From ea4917fda25cf1cfe86d8fc6dc40f80bd8f3cc46 Mon Sep 17 00:00:00 2001 From: Andrea Paris Date: Fri, 27 Feb 2026 20:06:03 +0100 Subject: [PATCH 08/12] gradient cleanup --- tsd/src/tsd/core/algorithms/vort.h | 125 +++++++++-------------------- 1 file changed, 39 insertions(+), 86 deletions(-) diff --git a/tsd/src/tsd/core/algorithms/vort.h b/tsd/src/tsd/core/algorithms/vort.h index 234432ad9..aa61bb7ba 100644 --- a/tsd/src/tsd/core/algorithms/vort.h +++ b/tsd/src/tsd/core/algorithms/vort.h @@ -53,90 +53,43 @@ static double q_crit(double j00, double j01, double j02, } // --------------------------------------------------------------------------- -// 1D finite-difference helpers — float* input, double* output, stride s. +// grad3D — finite differences along one axis of a 3D float field. // -// grad1D: uniform spacing; 2nd-order one-sided BCs, 2nd-order central. -// inv2h = 1/(2h). -// grad1D_nu: non-uniform spacing; 1st-order one-sided BCs, 2nd-order central. -// c[]: coordinate array (length n). -// --------------------------------------------------------------------------- - -static void grad1D( - const float *fc, double *gc, size_t n, size_t s, double inv2h) -{ - gc[0] = (-3.0*fc[0] + 4.0*fc[s] - fc[2*s]) * inv2h; - for (size_t i = 1; i < n-1; ++i) - gc[i*s] = ((double)fc[(i+1)*s] - (double)fc[(i-1)*s]) * inv2h; - gc[(n-1)*s] = ( 3.0*fc[(n-1)*s] - 4.0*fc[(n-2)*s] + fc[(n-3)*s]) * inv2h; -} - -static void grad1D_nu( - const float *fc, double *gc, size_t n, size_t s, const double *c) -{ - gc[0] = ((double)fc[s] - (double)fc[0]) / (c[1] - c[0]); - for (size_t i = 1; i < n-1; ++i) - gc[i*s] = ((double)fc[(i+1)*s] - (double)fc[(i-1)*s]) / (c[i+1] - c[i-1]); - gc[(n-1)*s] = ((double)fc[(n-1)*s] - (double)fc[(n-2)*s]) / (c[n-1] - c[n-2]); -} - -// --------------------------------------------------------------------------- -// Gradient helpers — float* input, double* output (scratch). -// All three take a coordinate array so uniform and rectilinear grids share the -// same interface; uniform fields pass the x_/y_/z_ arrays built from -// origin+spacing in extractStructuredRegular. +// axis 0 (x): stride 1, n=nx, outer=ny*nz passes +// axis 1 (y): stride nx, n=ny, outer=nx*nz passes +// axis 2 (z): stride nx*ny, n=nz, outer=nx*ny passes // -// Boundary accuracy: 1st-order one-sided. Interior: 2nd-order central. +// c[]: coordinate array along the axis (length n). +// Boundary: 1st-order one-sided. Interior: 2nd-order central. // --------------------------------------------------------------------------- - -// x is the fastest index, so every (z,y) pair is a contiguous row of nx -// floats — collapse both outer loops into a single row index. -inline void gradX(const float *u, - double *uGrad, - const double *x_, +static void grad3D(const float *fc, + double *gc, size_t nx, size_t ny, - size_t nz) -{ - const int nrows = (int)(ny * nz); - for (int r = 0; r < nrows; ++r) - grad1D_nu(u + (size_t)r*nx, uGrad + (size_t)r*nx, nx, 1, x_); -} - -// y has stride nx. Iterate over every (z,x) pair as the column base; -// the 1D y-pass then uses gc[y*nx] with no extra offset arithmetic. -inline void gradY(const float *v, - double *vGrad, - const double *y_, - size_t nx, - size_t ny, - size_t nz) -{ - const size_t slab = nx * ny; - for (int z = 0; z < (int)nz; ++z) - for (int x = 0; x < (int)nx; ++x) - grad1D_nu(v + (size_t)z*slab + x, vGrad + (size_t)z*slab + x, ny, nx, y_); -} - -// Non-uniform z spacing: 1st-order one-sided at boundaries, 2nd-order central -// in the interior. -// z has stride slab=nx*ny. Every (y,x) pair gives a z-column addressable as -// wc[k*slab] — flatten both outer loops into a single index over ny*nx. -inline void gradZ(const float *w, - double *wGrad, - const double *z_, - size_t nx, - size_t ny, - size_t nz) + size_t nz, + int axis, + const double *c) { const size_t slab = nx * ny; - const int ncols = (int)slab; - for (int yx = 0; yx < ncols; ++yx) { - const float *wc = w + yx; - double *gc = wGrad + yx; - gc[0] = ((double)wc[slab] - (double)wc[0]) / (z_[1] - z_[0]); - for (size_t k = 1; k < nz - 1; ++k) - gc[k * slab] = ((double)wc[(k+1)*slab] - (double)wc[(k-1)*slab]) / (z_[k+1] - z_[k-1]); - gc[(nz-1) * slab] = ((double)wc[(nz-1)*slab] - (double)wc[(nz-2)*slab]) / (z_[nz-1] - z_[nz-2]); + size_t n, s, outer; + if (axis == 0) { n = nx; s = 1; outer = ny * nz; } + else if (axis == 1) { n = ny; s = nx; outer = nx * nz; } + else { n = nz; s = slab; outer = slab; } + + auto pass = [n, s, c](const float *f, double *g) { + g[0] = ((double)f[s] - (double)f[0]) / (c[1] - c[0]); + for (size_t i = 1; i < n-1; ++i) + g[i*s] = ((double)f[(i+1)*s] - (double)f[(i-1)*s]) / (c[i+1] - c[i-1]); + g[(n-1)*s] = ((double)f[(n-1)*s] - (double)f[(n-2)*s]) / (c[n-1] - c[n-2]); + }; + + if (axis == 1) { + for (size_t z = 0; z < nz; ++z) + for (size_t x = 0; x < nx; ++x) + pass(fc + z*slab + x, gc + z*slab + x); + } else { + for (size_t r = 0; r < outer; ++r) + pass(fc + r*n*s, gc + r*n*s); } } @@ -232,17 +185,17 @@ inline void vort(const float *u, std::vector duy(len), dvy(len), dwy(len); std::vector duz(len), dvz(len), dwz(len); - gradX(u, dux.data(), x_, nx, ny, nz); - gradX(v, dvx.data(), x_, nx, ny, nz); - gradX(w, dwx.data(), x_, nx, ny, nz); + grad3D(u, dux.data(), nx, ny, nz, 0, x_); + grad3D(v, dvx.data(), nx, ny, nz, 0, x_); + grad3D(w, dwx.data(), nx, ny, nz, 0, x_); - gradY(u, duy.data(), y_, nx, ny, nz); - gradY(v, dvy.data(), y_, nx, ny, nz); - gradY(w, dwy.data(), y_, nx, ny, nz); + grad3D(u, duy.data(), nx, ny, nz, 1, y_); + grad3D(v, dvy.data(), nx, ny, nz, 1, y_); + grad3D(w, dwy.data(), nx, ny, nz, 1, y_); - gradZ(u, duz.data(), z_, nx, ny, nz); - gradZ(v, dvz.data(), z_, nx, ny, nz); - gradZ(w, dwz.data(), z_, nx, ny, nz); + grad3D(u, duz.data(), nx, ny, nz, 2, z_); + grad3D(v, dvz.data(), nx, ny, nz, 2, z_); + grad3D(w, dwz.data(), nx, ny, nz, 2, z_); vort_from_jacobians(u, v, From c7c7aa77826b8d90565d8c874ecb3f4288101825 Mon Sep 17 00:00:00 2001 From: Andrea Paris Date: Fri, 27 Feb 2026 20:18:32 +0100 Subject: [PATCH 09/12] computeVorticity cleanup --- .../tsd/io/procedural/computeVorticity.cpp | 167 +++++++----------- 1 file changed, 68 insertions(+), 99 deletions(-) diff --git a/tsd/src/tsd/io/procedural/computeVorticity.cpp b/tsd/src/tsd/io/procedural/computeVorticity.cpp index 09ad3d913..577a55166 100644 --- a/tsd/src/tsd/io/procedural/computeVorticity.cpp +++ b/tsd/src/tsd/io/procedural/computeVorticity.cpp @@ -47,6 +47,60 @@ struct FieldData size_t nx{0}, ny{0}, nz{0}; }; +// --------------------------------------------------------------------------- +// Shared helpers +// --------------------------------------------------------------------------- + +static void fillLinearCoords( + std::vector &dst, size_t n, double origin, double spacing) +{ + dst.resize(n); + for (size_t i = 0; i < n; ++i) + dst[i] = origin + i * spacing; +} + +static bool readCoordsParam(Scene &scene, + SpatialField *field, + const char *name, + std::vector &dst, + size_t n) +{ + auto *cp = field->parameter(name); + if (!cp || !anari::isArray(cp->value().type())) { + logError("[computeVorticity] field missing coord parameter '%s'", name); + return false; + } + auto ca = scene.getObject(cp->value().getAsObjectIndex()); + if (!ca) { + logError("[computeVorticity] field coord parameter '%s' is invalid", name); + return false; + } + dst.resize(n); + const float *src = ca->dataAs(); + for (size_t i = 0; i < n; ++i) + dst[i] = src[i]; + return true; +} + +static void rasterizeNanoVDB(const nanovdb::NanoGrid *grid, + nanovdb::Coord lo, + nanovdb::Coord hi, + FieldData &out) +{ + out.ownedData.resize(out.nx * out.ny * out.nz, 0.0f); + out.ptr = out.ownedData.data(); + auto acc = grid->getAccessor(); + for (int k = lo[2]; k <= hi[2]; ++k) { + for (int j = lo[1]; j <= hi[1]; ++j) { + float *row = out.ownedData.data() + + (size_t)(k - lo[2]) * out.ny * out.nx + + (size_t)(j - lo[1]) * out.nx; + for (int i = lo[0]; i <= hi[0]; ++i) + row[i - lo[0]] = acc.getValue(nanovdb::Coord(i, j, k)); + } + } +} + // --------------------------------------------------------------------------- // Extraction helpers // --------------------------------------------------------------------------- @@ -96,15 +150,9 @@ static bool extractStructuredRegular( // Point directly at the ANARI array's storage — no copy needed. out.ptr = arr->dataAs(); - out.x.resize(out.nx); - out.y.resize(out.ny); - out.z.resize(out.nz); - for (size_t i = 0; i < out.nx; ++i) - out.x[i] = origin.x + i * (double)spacing.x; - for (size_t j = 0; j < out.ny; ++j) - out.y[j] = origin.y + j * (double)spacing.y; - for (size_t k = 0; k < out.nz; ++k) - out.z[k] = origin.z + k * (double)spacing.z; + fillLinearCoords(out.x, out.nx, origin.x, spacing.x); + fillLinearCoords(out.y, out.ny, origin.y, spacing.y); + fillLinearCoords(out.z, out.nz, origin.z, spacing.z); return true; } @@ -156,30 +204,10 @@ static bool extractNanoVDB(Scene &scene, SpatialField *field, FieldData &out) double dy = (out.ny > 1) ? (worldHi[1] - worldLo[1]) / (out.ny - 1) : 1.0; double dz = (out.nz > 1) ? (worldHi[2] - worldLo[2]) / (out.nz - 1) : 1.0; - out.x.resize(out.nx); - out.y.resize(out.ny); - out.z.resize(out.nz); - for (size_t i = 0; i < out.nx; ++i) - out.x[i] = worldLo[0] + i * dx; - for (size_t j = 0; j < out.ny; ++j) - out.y[j] = worldLo[1] + j * dy; - for (size_t k = 0; k < out.nz; ++k) - out.z[k] = worldLo[2] + k * dz; - - // Rasterize sparse grid into a dense float buffer (no double conversion). - out.ownedData.resize(out.nx * out.ny * out.nz, 0.0f); - out.ptr = out.ownedData.data(); - auto acc = grid->getAccessor(); - - for (int k = lo[2]; k <= hi[2]; ++k) { - for (int j = lo[1]; j <= hi[1]; ++j) { - float *row = out.ownedData.data() - + (size_t)(k - lo[2]) * out.ny * out.nx - + (size_t)(j - lo[1]) * out.nx; - for (int i = lo[0]; i <= hi[0]; ++i) - row[i - lo[0]] = acc.getValue(nanovdb::Coord(i, j, k)); - } - } + fillLinearCoords(out.x, out.nx, worldLo[0], dx); + fillLinearCoords(out.y, out.ny, worldLo[1], dy); + fillLinearCoords(out.z, out.nz, worldLo[2], dz); + rasterizeNanoVDB(grid, lo, hi, out); return true; } @@ -220,32 +248,9 @@ static bool extractStructuredRectilinear( out.ptr = arr->dataAs(); - auto readCoords = [&](const char *name, - std::vector &dst, - size_t n) -> bool { - auto *cp = field->parameter(name); - if (!cp || !anari::isArray(cp->value().type())) { - logError( - "[computeVorticity] structuredRectilinear field missing '%s'", name); - return false; - } - auto ca = scene.getObject(cp->value().getAsObjectIndex()); - if (!ca) { - logError( - "[computeVorticity] structuredRectilinear '%s' array is invalid", - name); - return false; - } - dst.resize(n); - const float *src = ca->dataAs(); - for (size_t i = 0; i < n; ++i) - dst[i] = src[i]; - return true; - }; - - return readCoords("coordsX", out.x, out.nx) - && readCoords("coordsY", out.y, out.ny) - && readCoords("coordsZ", out.z, out.nz); + return readCoordsParam(scene, field, "coordsX", out.x, out.nx) + && readCoordsParam(scene, field, "coordsY", out.y, out.ny) + && readCoordsParam(scene, field, "coordsZ", out.z, out.nz); } static bool extractNanoVDBRectilinear( @@ -292,48 +297,12 @@ static bool extractNanoVDBRectilinear( return false; } - auto readCoords = [&](const char *name, - std::vector &dst, - size_t n) -> bool { - auto *cp = field->parameter(name); - if (!cp || !anari::isArray(cp->value().type())) { - logError( - "[computeVorticity] nanovdbRectilinear field missing '%s'", name); - return false; - } - auto ca = scene.getObject(cp->value().getAsObjectIndex()); - if (!ca) { - logError( - "[computeVorticity] nanovdbRectilinear '%s' array is invalid", name); - return false; - } - dst.resize(n); - const float *src = ca->dataAs(); - for (size_t i = 0; i < n; ++i) - dst[i] = src[i]; - return true; - }; - - if (!readCoords("coordsX", out.x, out.nx) - || !readCoords("coordsY", out.y, out.ny) - || !readCoords("coordsZ", out.z, out.nz)) + if (!readCoordsParam(scene, field, "coordsX", out.x, out.nx) + || !readCoordsParam(scene, field, "coordsY", out.y, out.ny) + || !readCoordsParam(scene, field, "coordsZ", out.z, out.nz)) return false; - // Rasterize sparse grid into a dense float buffer - out.ownedData.resize(out.nx * out.ny * out.nz, 0.0f); - out.ptr = out.ownedData.data(); - auto acc = grid->getAccessor(); - - for (int k = lo[2]; k <= hi[2]; ++k) { - for (int j = lo[1]; j <= hi[1]; ++j) { - float *row = out.ownedData.data() - + (size_t)(k - lo[2]) * out.ny * out.nx - + (size_t)(j - lo[1]) * out.nx; - for (int i = lo[0]; i <= hi[0]; ++i) - row[i - lo[0]] = acc.getValue(nanovdb::Coord(i, j, k)); - } - } - + rasterizeNanoVDB(grid, lo, hi, out); return true; } From 0faf222b97f73a0d899472f1e164cc8d52b0d6ee Mon Sep 17 00:00:00 2001 From: Andrea Paris Date: Fri, 27 Feb 2026 20:30:31 +0100 Subject: [PATCH 10/12] ui update, removed ... --- tsd/src/tsd/ui/imgui/Application.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tsd/src/tsd/ui/imgui/Application.cpp b/tsd/src/tsd/ui/imgui/Application.cpp index bb5809783..f6a625d65 100644 --- a/tsd/src/tsd/ui/imgui/Application.cpp +++ b/tsd/src/tsd/ui/imgui/Application.cpp @@ -341,7 +341,7 @@ void Application::uiMainMenuBar_Tools() ImGui::Separator(); - if (ImGui::MenuItem("Flow Analysis...")) + if (ImGui::MenuItem("Flow Analysis")) m_vorticityDialog->show(); ImGui::Separator(); From eada49661c3301b782e8aa51b989c15c76fe59c7 Mon Sep 17 00:00:00 2001 From: Andrea Paris Date: Fri, 27 Feb 2026 22:23:55 +0100 Subject: [PATCH 11/12] copyright update --- tsd/src/tsd/core/algorithms/vort.h | 2 +- tsd/src/tsd/ui/imgui/modals/VorticityDialog.cpp | 2 +- tsd/src/tsd/ui/imgui/modals/VorticityDialog.h | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tsd/src/tsd/core/algorithms/vort.h b/tsd/src/tsd/core/algorithms/vort.h index aa61bb7ba..1562faf2f 100644 --- a/tsd/src/tsd/core/algorithms/vort.h +++ b/tsd/src/tsd/core/algorithms/vort.h @@ -1,4 +1,4 @@ -// Copyright 2025-2026 NVIDIA Corporation +// Copyright 2026 NVIDIA Corporation // SPDX-License-Identifier: Apache-2.0 #pragma once diff --git a/tsd/src/tsd/ui/imgui/modals/VorticityDialog.cpp b/tsd/src/tsd/ui/imgui/modals/VorticityDialog.cpp index 8116e9be9..ddf132c57 100644 --- a/tsd/src/tsd/ui/imgui/modals/VorticityDialog.cpp +++ b/tsd/src/tsd/ui/imgui/modals/VorticityDialog.cpp @@ -1,4 +1,4 @@ -// Copyright 2025-2026 NVIDIA Corporation +// Copyright 2026 NVIDIA Corporation // SPDX-License-Identifier: Apache-2.0 #include "VorticityDialog.h" diff --git a/tsd/src/tsd/ui/imgui/modals/VorticityDialog.h b/tsd/src/tsd/ui/imgui/modals/VorticityDialog.h index 9351b734a..564ef1aa0 100644 --- a/tsd/src/tsd/ui/imgui/modals/VorticityDialog.h +++ b/tsd/src/tsd/ui/imgui/modals/VorticityDialog.h @@ -1,4 +1,4 @@ -// Copyright 2025-2026 NVIDIA Corporation +// Copyright 2026 NVIDIA Corporation // SPDX-License-Identifier: Apache-2.0 #pragma once From 26c45e65fe4511c1a235f5baa723b8a9bb7b378c Mon Sep 17 00:00:00 2001 From: Andrea Paris Date: Fri, 27 Feb 2026 22:33:25 +0100 Subject: [PATCH 12/12] formatter --- tsd/src/tsd/core/algorithms/vort.h | 141 ++++++++++----- tsd/src/tsd/io/importers/import_USD.cpp | 17 +- tsd/src/tsd/io/importers/import_VTI.cpp | 3 +- tsd/src/tsd/io/importers/import_VTU.cpp | 6 +- .../tsd/io/procedural/computeVorticity.cpp | 164 ++++++++++-------- .../tsd/ui/imgui/modals/VorticityDialog.cpp | 9 +- 6 files changed, 207 insertions(+), 133 deletions(-) diff --git a/tsd/src/tsd/core/algorithms/vort.h b/tsd/src/tsd/core/algorithms/vort.h index 1562faf2f..b08725cdd 100644 --- a/tsd/src/tsd/core/algorithms/vort.h +++ b/tsd/src/tsd/core/algorithms/vort.h @@ -10,46 +10,70 @@ // lambda2: middle eigenvalue of M = S^2 + O^2, where S=(J+J^T)/2, O=(J-J^T)/2. // M is symmetric with 6 unique entries computed as M = (J^2 + (J^T)^2) / 2. -// Eigenvalues via the trigonometric (Cardano) formula: eig_k = q + 2p*cos(phi + 2*pi*k/3). -// Middle eigenvalue recovered from the trace identity: eig1 = 3q - eig0 - eig2. -static double l2(double j00, double j01, double j02, - double j10, double j11, double j12, - double j20, double j21, double j22) +// Eigenvalues via the trigonometric (Cardano) formula: eig_k = q + 2p*cos(phi + +// 2*pi*k/3). Middle eigenvalue recovered from the trace identity: eig1 = 3q - +// eig0 - eig2. +static double l2(double j00, + double j01, + double j02, + double j10, + double j11, + double j12, + double j20, + double j21, + double j22) { // 6 unique entries of symmetric M = (J^2 + (J^T)^2) / 2 - const double m00 = j00*j00 + j01*j10 + j02*j20; - const double m11 = j10*j01 + j11*j11 + j12*j21; - const double m22 = j20*j02 + j21*j12 + j22*j22; - const double m01 = 0.5*(j00*j01 + j01*j11 + j02*j21 + j00*j10 + j10*j11 + j20*j12); - const double m02 = 0.5*(j00*j02 + j01*j12 + j02*j22 + j00*j20 + j10*j21 + j20*j22); - const double m12 = 0.5*(j10*j02 + j11*j12 + j12*j22 + j01*j20 + j11*j21 + j21*j22); - - const double q = (m00 + m11 + m22) / 3.0; - const double a = m00-q, d = m11-q, f = m22-q; - const double p1 = m01*m01 + m02*m02 + m12*m12; + const double m00 = j00 * j00 + j01 * j10 + j02 * j20; + const double m11 = j10 * j01 + j11 * j11 + j12 * j21; + const double m22 = j20 * j02 + j21 * j12 + j22 * j22; + const double m01 = 0.5 + * (j00 * j01 + j01 * j11 + j02 * j21 + j00 * j10 + j10 * j11 + j20 * j12); + const double m02 = 0.5 + * (j00 * j02 + j01 * j12 + j02 * j22 + j00 * j20 + j10 * j21 + j20 * j22); + const double m12 = 0.5 + * (j10 * j02 + j11 * j12 + j12 * j22 + j01 * j20 + j11 * j21 + j21 * j22); + + const double q = (m00 + m11 + m22) / 3.0; + const double a = m00 - q, d = m11 - q, f = m22 - q; + const double p1 = m01 * m01 + m02 * m02 + m12 * m12; if (p1 == 0.0) { // Diagonal — middle eigenvalue by inspection double e[3] = {a, d, f}; - if (e[0] < e[1]) std::swap(e[0], e[1]); - if (e[1] < e[2]) std::swap(e[1], e[2]); - if (e[0] < e[1]) std::swap(e[0], e[1]); + if (e[0] < e[1]) + std::swap(e[0], e[1]); + if (e[1] < e[2]) + std::swap(e[1], e[2]); + if (e[0] < e[1]) + std::swap(e[0], e[1]); return q + e[1]; } - const double p = std::sqrt((a*a + d*d + f*f + 2.0*p1) / 6.0); - const double r = (a*(d*f - m12*m12) - m01*(m01*f - m12*m02) + m02*(m01*m12 - d*m02)) - / (2.0*p*p*p); + const double p = std::sqrt((a * a + d * d + f * f + 2.0 * p1) / 6.0); + const double r = (a * (d * f - m12 * m12) - m01 * (m01 * f - m12 * m02) + + m02 * (m01 * m12 - d * m02)) + / (2.0 * p * p * p); const double phi = std::acos(std::max(-1.0, std::min(1.0, r))) / 3.0; - const double e0 = q + 2.0*p*std::cos(phi); - const double e2 = q + 2.0*p*std::cos(phi + 2.0943951023931953); // phi + 2*pi/3 - return 3.0*q - e0 - e2; // middle eigenvalue via trace identity + const double e0 = q + 2.0 * p * std::cos(phi); + const double e2 = + q + 2.0 * p * std::cos(phi + 2.0943951023931953); // phi + 2*pi/3 + return 3.0 * q - e0 - e2; // middle eigenvalue via trace identity } -// Q-criterion: 0.5*(||O||^2 - ||S||^2) = -0.5*tr(J^2) = -0.5*(j00^2+j11^2+j22^2 + 2*(j01*j10+j02*j20+j12*j21)) -static double q_crit(double j00, double j01, double j02, - double j10, double j11, double j12, - double j20, double j21, double j22) +// Q-criterion: 0.5*(||O||^2 - ||S||^2) = -0.5*tr(J^2) = -0.5*(j00^2+j11^2+j22^2 +// + 2*(j01*j10+j02*j20+j12*j21)) +static double q_crit(double j00, + double j01, + double j02, + double j10, + double j11, + double j12, + double j20, + double j21, + double j22) { - return -0.5*(j00*j00 + j11*j11 + j22*j22 + 2.0*(j01*j10 + j02*j20 + j12*j21)); + return -0.5 + * (j00 * j00 + j11 * j11 + j22 * j22 + + 2.0 * (j01 * j10 + j02 * j20 + j12 * j21)); } // --------------------------------------------------------------------------- @@ -72,24 +96,36 @@ static void grad3D(const float *fc, { const size_t slab = nx * ny; size_t n, s, outer; - if (axis == 0) { n = nx; s = 1; outer = ny * nz; } - else if (axis == 1) { n = ny; s = nx; outer = nx * nz; } - else { n = nz; s = slab; outer = slab; } + if (axis == 0) { + n = nx; + s = 1; + outer = ny * nz; + } else if (axis == 1) { + n = ny; + s = nx; + outer = nx * nz; + } else { + n = nz; + s = slab; + outer = slab; + } auto pass = [n, s, c](const float *f, double *g) { - g[0] = ((double)f[s] - (double)f[0]) / (c[1] - c[0]); - for (size_t i = 1; i < n-1; ++i) - g[i*s] = ((double)f[(i+1)*s] - (double)f[(i-1)*s]) / (c[i+1] - c[i-1]); - g[(n-1)*s] = ((double)f[(n-1)*s] - (double)f[(n-2)*s]) / (c[n-1] - c[n-2]); + g[0] = ((double)f[s] - (double)f[0]) / (c[1] - c[0]); + for (size_t i = 1; i < n - 1; ++i) + g[i * s] = ((double)f[(i + 1) * s] - (double)f[(i - 1) * s]) + / (c[i + 1] - c[i - 1]); + g[(n - 1) * s] = ((double)f[(n - 1) * s] - (double)f[(n - 2) * s]) + / (c[n - 1] - c[n - 2]); }; if (axis == 1) { for (size_t z = 0; z < nz; ++z) for (size_t x = 0; x < nx; ++x) - pass(fc + z*slab + x, gc + z*slab + x); + pass(fc + z * slab + x, gc + z * slab + x); } else { for (size_t r = 0; r < outer; ++r) - pass(fc + r*n*s, gc + r*n*s); + pass(fc + r * n * s, gc + r * n * s); } } @@ -98,7 +134,8 @@ static void grad3D(const float *fc, // components. // // Inputs: u, v, w — float velocity components (any layout, length len) -// dux..dwz — Jacobian entries: d(u,v,w)/d(x,y,z), double, length len +// dux..dwz — Jacobian entries: d(u,v,w)/d(x,y,z), double, length +// len // naming: d, e.g. duy = ∂u/∂y // Outputs: vorticity, helicity, lambda2, qCriterion — float*, written in-place. // Any output pointer may be null; null outputs are simply skipped. @@ -140,13 +177,27 @@ inline void vort_from_jacobians(const float *u, } if (lambda2 || qCriterion) { if (lambda2) - lambda2[i] = (float)(-std::min(l2(dux[i], duy[i], duz[i], - dvx[i], dvy[i], dvz[i], - dwx[i], dwy[i], dwz[i]), 0.0)); + lambda2[i] = (float)(-std::min(l2(dux[i], + duy[i], + duz[i], + dvx[i], + dvy[i], + dvz[i], + dwx[i], + dwy[i], + dwz[i]), + 0.0)); if (qCriterion) - qCriterion[i] = (float)std::max(q_crit(dux[i], duy[i], duz[i], - dvx[i], dvy[i], dvz[i], - dwx[i], dwy[i], dwz[i]), 0.0); + qCriterion[i] = (float)std::max(q_crit(dux[i], + duy[i], + duz[i], + dvx[i], + dvy[i], + dvz[i], + dwx[i], + dwy[i], + dwz[i]), + 0.0); } } } diff --git a/tsd/src/tsd/io/importers/import_USD.cpp b/tsd/src/tsd/io/importers/import_USD.cpp index f124fc868..782cea05a 100644 --- a/tsd/src/tsd/io/importers/import_USD.cpp +++ b/tsd/src/tsd/io/importers/import_USD.cpp @@ -21,9 +21,9 @@ #include #include #include +#include #include #include -#include #include #include #include @@ -679,8 +679,8 @@ static void import_usd_points(Scene &scene, pointsPrim.GetPointsAttr().GetTimeSamples(&timeSamples); // Build position+radius arrays for one time step - auto buildFrame = - [&](pxr::UsdTimeCode tc) -> std::pair, ObjectUsePtr> { + auto buildFrame = [&](pxr::UsdTimeCode tc) + -> std::pair, ObjectUsePtr> { pxr::VtArray pts; pxr::VtArray wids; pointsPrim.GetPointsAttr().Get(&pts, tc); @@ -1471,9 +1471,7 @@ static void import_usd_prim_recursive(Scene &scene, } } -void import_USD(Scene &scene, - const char *filepath, - LayerNodeRef location) +void import_USD(Scene &scene, const char *filepath, LayerNodeRef location) { pxr::UsdStageRefPtr stage = pxr::UsdStage::Open(filepath); if (!stage) { @@ -1504,14 +1502,13 @@ void import_USD(Scene &scene, for (pxr::UsdPrim const &prim : stage->Traverse()) { // if (prim.IsPrototype()) continue; if (prim.GetParent() && prim.GetParent().IsPseudoRoot()) { - import_usd_prim_recursive(scene, prim, usd_root, xformCache, basePath, pxr::GfMatrix4d(1.0)); + import_usd_prim_recursive( + scene, prim, usd_root, xformCache, basePath, pxr::GfMatrix4d(1.0)); } } } #else -void import_USD(Scene &scene, - const char *filepath, - LayerNodeRef location) +void import_USD(Scene &scene, const char *filepath, LayerNodeRef location) { tsd::core::logError("[import_USD] USD not enabled in TSD build."); } diff --git a/tsd/src/tsd/io/importers/import_VTI.cpp b/tsd/src/tsd/io/importers/import_VTI.cpp index e0d9b1f89..38306505d 100644 --- a/tsd/src/tsd/io/importers/import_VTI.cpp +++ b/tsd/src/tsd/io/importers/import_VTI.cpp @@ -74,8 +74,7 @@ SpatialFieldRef import_VTI(Scene &scene, const char *filepath) std::string nameY = baseName + "_y"; std::string nameZ = baseName + "_z"; - size_t n = - (size_t)dims[0] * (size_t)dims[1] * (size_t)dims[2]; + size_t n = (size_t)dims[0] * (size_t)dims[1] * (size_t)dims[2]; auto arrX = scene.createArray(ANARI_FLOAT32, dims[0], dims[1], dims[2]); auto arrY = scene.createArray(ANARI_FLOAT32, dims[0], dims[1], dims[2]); diff --git a/tsd/src/tsd/io/importers/import_VTU.cpp b/tsd/src/tsd/io/importers/import_VTU.cpp index 40de88214..58a214760 100644 --- a/tsd/src/tsd/io/importers/import_VTU.cpp +++ b/tsd/src/tsd/io/importers/import_VTU.cpp @@ -351,9 +351,9 @@ static SpatialFieldRef createFieldFromVolumeCells( if (!array) continue; int nComp = array->GetNumberOfComponents(); - std::string arrName = - (array->GetName() && array->GetName()[0] != '\0') ? array->GetName() - : baseName; + std::string arrName = (array->GetName() && array->GetName()[0] != '\0') + ? array->GetName() + : baseName; if (nComp == 1) { auto dataArr = scene.createArray(ANARI_FLOAT32, numPoints); diff --git a/tsd/src/tsd/io/procedural/computeVorticity.cpp b/tsd/src/tsd/io/procedural/computeVorticity.cpp index 577a55166..8d39f9ce0 100644 --- a/tsd/src/tsd/io/procedural/computeVorticity.cpp +++ b/tsd/src/tsd/io/procedural/computeVorticity.cpp @@ -41,7 +41,8 @@ using namespace tsd::core; struct FieldData { - const float *ptr{nullptr}; // non-owning (structuredRegular) or alias into ownedData + const float *ptr{ + nullptr}; // non-owning (structuredRegular) or alias into ownedData std::vector ownedData; // owns data for NanoVDB rasterization std::vector x, y, z; // world-space coordinate arrays size_t nx{0}, ny{0}, nz{0}; @@ -92,8 +93,7 @@ static void rasterizeNanoVDB(const nanovdb::NanoGrid *grid, auto acc = grid->getAccessor(); for (int k = lo[2]; k <= hi[2]; ++k) { for (int j = lo[1]; j <= hi[1]; ++j) { - float *row = out.ownedData.data() - + (size_t)(k - lo[2]) * out.ny * out.nx + float *row = out.ownedData.data() + (size_t)(k - lo[2]) * out.ny * out.nx + (size_t)(j - lo[1]) * out.nx; for (int i = lo[0]; i <= hi[0]; ++i) row[i - lo[0]] = acc.getValue(nanovdb::Coord(i, j, k)); @@ -140,7 +140,8 @@ static bool extractStructuredRegular( auto *p_orig = field->parameter("origin"); auto *p_spc = field->parameter("spacing"); if (!p_orig || !p_spc) { - logError("[computeVorticity] structuredRegular field missing origin/spacing"); + logError( + "[computeVorticity] structuredRegular field missing origin/spacing"); return false; } @@ -172,8 +173,7 @@ static bool extractNanoVDB(Scene &scene, SpatialField *field, FieldData &out) } const uint8_t *rawData = static_cast(arr->data()); - const auto *meta = - reinterpret_cast(rawData); + const auto *meta = reinterpret_cast(rawData); if (meta->gridType() != nanovdb::GridType::Float) { logError( @@ -258,21 +258,18 @@ static bool extractNanoVDBRectilinear( { auto *p = field->parameter("data"); if (!p || !anari::isArray(p->value().type())) { - logError( - "[computeVorticity] nanovdbRectilinear field has no 'data' array"); + logError("[computeVorticity] nanovdbRectilinear field has no 'data' array"); return false; } auto arr = scene.getObject(p->value().getAsObjectIndex()); if (!arr || !arr->data()) { - logError( - "[computeVorticity] nanovdbRectilinear 'data' array is invalid"); + logError("[computeVorticity] nanovdbRectilinear 'data' array is invalid"); return false; } const uint8_t *rawData = static_cast(arr->data()); - const auto *meta = - reinterpret_cast(rawData); + const auto *meta = reinterpret_cast(rawData); if (meta->gridType() != nanovdb::GridType::Float) { logError( @@ -347,8 +344,8 @@ static VolumeRef wrapAsVolume(Scene &scene, const math::float3 &spacing, LayerNodeRef location) { - auto field = - scene.createObject(tokens::spatial_field::structuredRegular); + auto field = scene.createObject( + tokens::spatial_field::structuredRegular); field->setName(name.c_str()); field->setParameter("origin", origin); field->setParameter("spacing", spacing); @@ -378,7 +375,6 @@ static VolumeRef wrapAsVolume(Scene &scene, #if TSD_USE_VTK - static VorticityResult computeVorticityUnstructured(Scene &scene, SpatialField *u, SpatialField *v, @@ -391,8 +387,7 @@ static VorticityResult computeVorticityUnstructured(Scene &scene, // Step 1: numPoints from u's vertex.position array auto *posParam = u->parameter("vertex.position"); if (!posParam || !anari::isArray(posParam->value().type())) { - logError( - "[computeVorticity] unstructured field missing 'vertex.position'"); + logError("[computeVorticity] unstructured field missing 'vertex.position'"); return result; } auto posArr = scene.getObject(posParam->value().getAsObjectIndex()); @@ -484,8 +479,9 @@ static VorticityResult computeVorticityUnstructured(Scene &scene, std::vector ptIds(npts); for (size_t j = 0; j < npts; ++j) ptIds[j] = static_cast(connectivity[start + j]); - vgrid->InsertNextCell( - static_cast(cellTypeData[ci]), static_cast(npts), ptIds.data()); + vgrid->InsertNextCell(static_cast(cellTypeData[ci]), + static_cast(npts), + ptIds.data()); } // Add interleaved velocity array (u,v,w per point) @@ -508,17 +504,14 @@ static VorticityResult computeVorticityUnstructured(Scene &scene, gradFilter->SetComputeGradient(1); gradFilter->Update(); - auto *gradOutput = - vtkUnstructuredGrid::SafeDownCast(gradFilter->GetOutput()); + auto *gradOutput = vtkUnstructuredGrid::SafeDownCast(gradFilter->GetOutput()); if (!gradOutput) { logError("[computeVorticity] vtkGradientFilter failed"); return result; } - vtkDataArray *gradData = - gradOutput->GetPointData()->GetArray("gradients"); + vtkDataArray *gradData = gradOutput->GetPointData()->GetArray("gradients"); if (!gradData || gradData->GetNumberOfComponents() != 9) { - logError( - "[computeVorticity] gradient filter output has unexpected format"); + logError("[computeVorticity] gradient filter output has unexpected format"); return result; } @@ -533,9 +526,15 @@ static VorticityResult computeVorticityUnstructured(Scene &scene, for (size_t i = 0; i < numPoints; ++i) { double t[9]; gradData->GetTuple(static_cast(i), t); - dux[i] = t[0]; duy[i] = t[1]; duz[i] = t[2]; - dvx[i] = t[3]; dvy[i] = t[4]; dvz[i] = t[5]; - dwx[i] = t[6]; dwy[i] = t[7]; dwz[i] = t[8]; + dux[i] = t[0]; + duy[i] = t[1]; + duz[i] = t[2]; + dvx[i] = t[3]; + dvy[i] = t[4]; + dvz[i] = t[5]; + dwx[i] = t[6]; + dwy[i] = t[7]; + dwz[i] = t[8]; } // Step 7: allocate output arrays and compute @@ -545,25 +544,38 @@ static VorticityResult computeVorticityUnstructured(Scene &scene, auto arr = scene.createArray(ANARI_FLOAT32, numPoints); return {arr, arr->mapAs()}; }; - auto [lambda2Arr, lambda2Out] = makeOutBuf(opts.lambda2); - auto [qCritArr, qCritOut] = makeOutBuf(opts.qCriterion); + auto [lambda2Arr, lambda2Out] = makeOutBuf(opts.lambda2); + auto [qCritArr, qCritOut] = makeOutBuf(opts.qCriterion); auto [vorticityArr, vorticityOut] = makeOutBuf(opts.vorticity); - auto [helicityArr, helicityOut] = makeOutBuf(opts.helicity); + auto [helicityArr, helicityOut] = makeOutBuf(opts.helicity); vort_from_jacobians(uPtr, vPtr, wPtr, - dux.data(), dvx.data(), dwx.data(), - duy.data(), dvy.data(), dwy.data(), - duz.data(), dvz.data(), dwz.data(), - vorticityOut, helicityOut, lambda2Out, qCritOut, + dux.data(), + dvx.data(), + dwx.data(), + duy.data(), + dvy.data(), + dwy.data(), + duz.data(), + dvz.data(), + dwz.data(), + vorticityOut, + helicityOut, + lambda2Out, + qCritOut, numPoints); // Step 8: unmap - if (lambda2Arr) lambda2Arr->unmap(); - if (qCritArr) qCritArr->unmap(); - if (vorticityArr) vorticityArr->unmap(); - if (helicityArr) helicityArr->unmap(); + if (lambda2Arr) + lambda2Arr->unmap(); + if (qCritArr) + qCritArr->unmap(); + if (vorticityArr) + vorticityArr->unmap(); + if (helicityArr) + helicityArr->unmap(); // Step 9: Resample onto a structured regular grid and create volumes. // VisRTX does not support the "unstructured" SpatialField subtype, so we @@ -584,10 +596,10 @@ static VorticityResult computeVorticityUnstructured(Scene &scene, vtkarr->SetValue(static_cast(i), src[i]); vgrid->GetPointData()->AddArray(vtkarr); }; - addToGrid("lambda2", lambda2Arr, opts.lambda2); - addToGrid("q_criterion", qCritArr, opts.qCriterion); - addToGrid("vorticity", vorticityArr, opts.vorticity); - addToGrid("helicity", helicityArr, opts.helicity); + addToGrid("lambda2", lambda2Arr, opts.lambda2); + addToGrid("q_criterion", qCritArr, opts.qCriterion); + addToGrid("vorticity", vorticityArr, opts.vorticity); + addToGrid("helicity", helicityArr, opts.helicity); // 9b. Compute output grid resolution proportional to bounding box extents. double bounds[6]; @@ -601,20 +613,23 @@ static VorticityResult computeVorticityUnstructured(Scene &scene, return result; } const int MAX_RES = 64; - const size_t resX = std::max(size_t(2), size_t(std::round(MAX_RES * bx / maxB))); - const size_t resY = std::max(size_t(2), size_t(std::round(MAX_RES * by / maxB))); - const size_t resZ = std::max(size_t(2), size_t(std::round(MAX_RES * bz / maxB))); + const size_t resX = + std::max(size_t(2), size_t(std::round(MAX_RES * bx / maxB))); + const size_t resY = + std::max(size_t(2), size_t(std::round(MAX_RES * by / maxB))); + const size_t resZ = + std::max(size_t(2), size_t(std::round(MAX_RES * bz / maxB))); logStatus("[computeVorticity] resampling to %zu×%zu×%zu structured grid...", - resX, resY, resZ); + resX, + resY, + resZ); auto imData = vtkSmartPointer::New(); - imData->SetDimensions(static_cast(resX), - static_cast(resY), - static_cast(resZ)); + imData->SetDimensions( + static_cast(resX), static_cast(resY), static_cast(resZ)); imData->SetOrigin(bounds[0], bounds[2], bounds[4]); - imData->SetSpacing( - bx / (resX - 1), by / (resY - 1), bz / (resZ - 1)); + imData->SetSpacing(bx / (resX - 1), by / (resY - 1), bz / (resZ - 1)); auto probe = vtkSmartPointer::New(); probe->SetSourceData(vgrid); @@ -631,8 +646,7 @@ static VorticityResult computeVorticityUnstructured(Scene &scene, const size_t resTotal = resX * resY * resZ; const math::float3 origin{ (float)bounds[0], (float)bounds[2], (float)bounds[4]}; - const math::float3 spacing{ - (float)(bx / (resX - 1)), + const math::float3 spacing{(float)(bx / (resX - 1)), (float)(by / (resY - 1)), (float)(bz / (resZ - 1))}; @@ -652,10 +666,10 @@ static VorticityResult computeVorticityUnstructured(Scene &scene, }; logStatus("[computeVorticity] creating output volumes..."); - result.lambda2 = wrapResampled("lambda2", opts.lambda2); + result.lambda2 = wrapResampled("lambda2", opts.lambda2); result.qCriterion = wrapResampled("q_criterion", opts.qCriterion); - result.vorticity = wrapResampled("vorticity", opts.vorticity); - result.helicity = wrapResampled("helicity", opts.helicity); + result.vorticity = wrapResampled("vorticity", opts.vorticity); + result.helicity = wrapResampled("helicity", opts.helicity); logStatus("[computeVorticity] done."); return result; @@ -692,14 +706,19 @@ VorticityResult computeVorticity(Scene &scene, return result; if (uData.nx != vData.nx || uData.nx != wData.nx || uData.ny != vData.ny - || uData.ny != wData.ny || uData.nz != vData.nz - || uData.nz != wData.nz) { + || uData.ny != wData.ny || uData.nz != vData.nz || uData.nz != wData.nz) { logError( "[computeVorticity] U/V/W fields have mismatched dimensions: " "U=(%zu,%zu,%zu) V=(%zu,%zu,%zu) W=(%zu,%zu,%zu)", - uData.nx, uData.ny, uData.nz, - vData.nx, vData.ny, vData.nz, - wData.nx, wData.ny, wData.nz); + uData.nx, + uData.ny, + uData.nz, + vData.nx, + vData.ny, + vData.nz, + wData.nx, + wData.ny, + wData.nz); return result; } @@ -707,7 +726,9 @@ VorticityResult computeVorticity(Scene &scene, logStatus( "[computeVorticity] computing vortical quantities on %zux%zux%zu grid...", - nx, ny, nz); + nx, + ny, + nz); // Pre-allocate one ANARI float32 array per selected output and map it so // vort() can write directly into the final storage — no intermediate buffers. @@ -717,10 +738,10 @@ VorticityResult computeVorticity(Scene &scene, auto arr = scene.createArray(ANARI_FLOAT32, nx, ny, nz); return {arr, arr->mapAs()}; }; - auto [lambda2Arr, lambda2Out] = makeOutBuf(opts.lambda2); - auto [qCritArr, qCritOut] = makeOutBuf(opts.qCriterion); + auto [lambda2Arr, lambda2Out] = makeOutBuf(opts.lambda2); + auto [qCritArr, qCritOut] = makeOutBuf(opts.qCriterion); auto [vorticityArr, vorticityOut] = makeOutBuf(opts.vorticity); - auto [helicityArr, helicityOut] = makeOutBuf(opts.helicity); + auto [helicityArr, helicityOut] = makeOutBuf(opts.helicity); // Compute — vort() reads float* velocity directly, writes float* outputs. // Null output pointers are silently skipped inside vort(). @@ -763,8 +784,15 @@ VorticityResult computeVorticity(Scene &scene, result.qCriterion = wrapAsVolume( scene, "q_criterion", qCritArr, nx, ny, nz, origin, spacing, location); if (vorticityArr) - result.vorticity = wrapAsVolume( - scene, "vorticity", vorticityArr, nx, ny, nz, origin, spacing, location); + result.vorticity = wrapAsVolume(scene, + "vorticity", + vorticityArr, + nx, + ny, + nz, + origin, + spacing, + location); if (helicityArr) result.helicity = wrapAsVolume( scene, "helicity", helicityArr, nx, ny, nz, origin, spacing, location); diff --git a/tsd/src/tsd/ui/imgui/modals/VorticityDialog.cpp b/tsd/src/tsd/ui/imgui/modals/VorticityDialog.cpp index ddf132c57..d99ccfcae 100644 --- a/tsd/src/tsd/ui/imgui/modals/VorticityDialog.cpp +++ b/tsd/src/tsd/ui/imgui/modals/VorticityDialog.cpp @@ -12,13 +12,12 @@ // imgui #include "imgui.h" // std -#include #include +#include namespace tsd::ui::imgui { -VorticityDialog::VorticityDialog(Application *app) - : Modal(app, "Flow Analysis") +VorticityDialog::VorticityDialog(Application *app) : Modal(app, "Flow Analysis") {} void VorticityDialog::buildUI() @@ -96,8 +95,8 @@ void VorticityDialog::buildUI() ImGui::SameLine(); - bool canRun = (m_uIdx >= 0 && m_vIdx >= 0 && m_wIdx >= 0 - && m_uIdx != m_vIdx && m_uIdx != m_wIdx && m_vIdx != m_wIdx); + bool canRun = (m_uIdx >= 0 && m_vIdx >= 0 && m_wIdx >= 0 && m_uIdx != m_vIdx + && m_uIdx != m_wIdx && m_vIdx != m_wIdx); ImGui::BeginDisabled(!canRun);