diff --git a/tsd/src/tsd/core/algorithms/vort.h b/tsd/src/tsd/core/algorithms/vort.h new file mode 100644 index 000000000..b08725cdd --- /dev/null +++ b/tsd/src/tsd/core/algorithms/vort.h @@ -0,0 +1,270 @@ +// Copyright 2026 NVIDIA Corporation +// SPDX-License-Identifier: Apache-2.0 + +#pragma once + +#include +#include +#include +#include + +// 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) +{ + // 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; + 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 - 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||^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)); +} + +// --------------------------------------------------------------------------- +// grad3D — finite differences along one axis of a 3D float field. +// +// 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 +// +// c[]: coordinate array along the axis (length n). +// Boundary: 1st-order one-sided. Interior: 2nd-order central. +// --------------------------------------------------------------------------- +static void grad3D(const float *fc, + double *gc, + size_t nx, + size_t ny, + size_t nz, + int axis, + const double *c) +{ + 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; + } + + 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); + } +} + +// --------------------------------------------------------------------------- +// 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) { + 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)); + 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); + } + } +} + +// --------------------------------------------------------------------------- +// 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. +// +// 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) +{ + if (!vorticity && !helicity && !lambda2 && !qCriterion) + return; + + 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); + + 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_); + + 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_); + + 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, + 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/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_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 1270238ca..38306505d 100644 --- a/tsd/src/tsd/io/importers/import_VTI.cpp +++ b/tsd/src/tsd/io/importers/import_VTI.cpp @@ -51,25 +51,87 @@ 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/importers/import_VTU.cpp b/tsd/src/tsd/io/importers/import_VTU.cpp index 50e90e8fb..58a214760 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 @@ -411,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.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..8d39f9ce0 --- /dev/null +++ b/tsd/src/tsd/io/procedural/computeVorticity.cpp @@ -0,0 +1,804 @@ +// 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 +#if TSD_USE_VTK +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#endif + +namespace tsd::io { + +using namespace tsd::core; + +// --------------------------------------------------------------------------- +// 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 +{ + 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}; +}; + +// --------------------------------------------------------------------------- +// 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 +// --------------------------------------------------------------------------- + +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(); + + // Point directly at the ANARI array's storage — no copy needed. + out.ptr = arr->dataAs(); + + 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; +} + +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; + } + + 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; + + 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; +} + +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(); + + 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( + 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; + } + + 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; + + rasterizeNanoVDB(grid, lo, hi, out); + 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::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, structuredRectilinear, nanovdb, and " + "nanovdbRectilinear are supported (unstructured and amr require " + "mesh interpolation first)", + field->subtype().c_str()); + return false; + } +} + +// --------------------------------------------------------------------------- +// 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, + ArrayRef &dataArr, + size_t nx, + size_t ny, + size_t nz, + const math::float3 &origin, + const math::float3 &spacing, + LayerNodeRef location) +{ + auto field = scene.createObject( + tokens::spatial_field::structuredRegular); + field->setName(name.c_str()); + field->setParameter("origin", origin); + field->setParameter("spacing", spacing); + field->setParameterObject("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; +} + +// --------------------------------------------------------------------------- +// Unstructured path — only compiled when VTK is available +// --------------------------------------------------------------------------- + +#if TSD_USE_VTK + +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 auto *rawPos = posArr->dataAs(); + for (size_t i = 0; i < numPoints; ++i) + points->SetPoint( + static_cast(i), rawPos[i].x, rawPos[i].y, rawPos[i].z); + 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(); + + // 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; + } + + // 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; +} + +#endif // TSD_USE_VTK + +// --------------------------------------------------------------------------- +// Public API +// --------------------------------------------------------------------------- + +VorticityResult computeVorticity(Scene &scene, + SpatialField *u, + SpatialField *v, + SpatialField *w, + LayerNodeRef location, + VorticityOptions opts) +{ + 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; + if (!extractFieldData(scene, u, uData)) + return result; + if (!extractFieldData(scene, v, vData)) + return result; + if (!extractFieldData(scene, w, wData)) + 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) { + 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; + } + + 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); + + // Pre-allocate one ANARI float32 array per selected output and map it so + // vort() can write directly into the final storage — no intermediate buffers. + 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(). + 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..."); + + 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 (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; +} + +} // 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..f6a625d65 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..d99ccfcae --- /dev/null +++ b/tsd/src/tsd/ui/imgui/modals/VorticityDialog.cpp @@ -0,0 +1,134 @@ +// Copyright 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..564ef1aa0 --- /dev/null +++ b/tsd/src/tsd/ui/imgui/modals/VorticityDialog.h @@ -0,0 +1,31 @@ +// Copyright 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