Skip to content

Commit 4f477be

Browse files
tarcilatarcila
committed
SILO importer to use ROI for ghost zones
1 parent ebb1427 commit 4f477be

File tree

1 file changed

+97
-130
lines changed

1 file changed

+97
-130
lines changed

‎tsd/src/tsd/io/importers/import_SILO.cpp‎

Lines changed: 97 additions & 130 deletions
Original file line numberDiff line numberDiff line change
@@ -7,9 +7,13 @@
77

88
#include "tsd/core/ColorMapUtil.hpp"
99
#include "tsd/core/Logging.hpp"
10+
#include "tsd/core/TSDMath.hpp"
1011
#include "tsd/io/importers.hpp"
1112
#include "tsd/io/importers/detail/importer_common.hpp"
1213

14+
// anari
15+
#include <anari/anari_cpp/Traits.h>
16+
1317
#if TSD_USE_SILO
1418
// silo
1519
#include <silo.h>
@@ -23,6 +27,10 @@
2327
#endif
2428
#endif
2529

30+
namespace anari {
31+
ANARI_TYPEFOR_SPECIALIZATION(tsd::math::box3, ANARI_FLOAT32_BOX3);
32+
} // namespace anari
33+
2634
namespace tsd::io {
2735

2836
using namespace tsd::core;
@@ -419,162 +427,108 @@ static SpatialFieldRef createFieldFromQuadMesh(Scene &scene,
419427
tokens::spatial_field::structuredRegular);
420428
field->setName(fieldName.c_str());
421429

422-
// Check for ghost zones using ghost_zone_labels
423-
// Note: min_index/max_index refer to NODE indices, not zone indices!
424-
int3 realNodeFirst(0, 0, 0);
425-
int3 realNodeLast(mesh->dims[0] - 1, mesh->dims[1] - 1, mesh->dims[2] - 1);
426-
bool hasGhostZones = false;
427-
428-
if (mesh->ghost_zone_labels) {
429-
hasGhostZones = true;
430-
logStatus("[import_SILO] mesh '%s' has ghost zone labels", meshName);
431-
432-
// Use min_index/max_index to determine real NODE bounds
433-
// These define the extent of real (non-ghost) nodes
434-
for (int i = 0; i < mesh->ndims && i < 3; i++) {
435-
realNodeFirst[i] = mesh->min_index[i];
436-
realNodeLast[i] = mesh->max_index[i];
437-
}
438-
logStatus("[import_SILO] real nodes: [%d:%d, %d:%d, %d:%d]",
439-
realNodeFirst.x,
440-
realNodeLast.x,
441-
realNodeFirst.y,
442-
realNodeLast.y,
443-
realNodeFirst.z,
444-
realNodeLast.z);
445-
}
446-
447-
// Get dimensions (only real zones)
448430
// Zones are defined between nodes, so zone count = node_count - 1
449-
int3 totalNodeDims(mesh->dims[0], mesh->dims[1], mesh->dims[2]);
450-
int3 totalDims(mesh->dims[0] - 1, mesh->dims[1] - 1, mesh->dims[2] - 1);
451-
int3 realZoneFirst(0, 0, 0);
452-
int3 realZoneLast = totalDims;
453-
int3 dims;
431+
int3 nodes(mesh->dims[0], mesh->dims[1], mesh->dims[2]);
432+
int3 zones(mesh->dims[0] - 1, mesh->dims[1] - 1, mesh->dims[2] - 1);
454433

455-
if (hasGhostZones) {
456-
// Real zones are between real nodes
457-
// If real nodes go from index A to B, real zones go from A to B-1
458-
realZoneFirst = realNodeFirst;
459-
realZoneLast = int3(std::max(realNodeLast.x - 1, realNodeFirst.x),
460-
std::max(realNodeLast.y - 1, realNodeFirst.y),
461-
std::max(realNodeLast.z - 1, realNodeFirst.z));
462-
463-
dims = int3(realZoneLast.x - realZoneFirst.x + 1,
464-
realZoneLast.y - realZoneFirst.y + 1,
465-
realZoneLast.z - realZoneFirst.z + 1);
434+
logStatus(
435+
"[import_SILO] total node dims:"
436+
" %d x %d x %d, total zone dims: %d x %d x %d",
437+
nodes.x,
438+
nodes.y,
439+
nodes.z,
440+
zones.x,
441+
zones.y,
442+
zones.z);
443+
444+
// Check for ghost zones, keeping only phony cells.
445+
// Phony is a term coming from the SILO headers.
446+
int3 minIndex(mesh->min_index[0], mesh->min_index[1], mesh->min_index[2]);
447+
int3 maxIndex(mesh->max_index[0], mesh->max_index[1], mesh->max_index[2]);
448+
bool hasGhostZones = minIndex != int3(0)
449+
|| maxIndex
450+
!= int3(mesh->dims[0] - 1, mesh->dims[1] - 1, mesh->dims[2] - 1);
451+
452+
// Check if data is node-centered or cell-centered
453+
bool isNodeCentered = (var->centering == DB_NODECENT);
466454

467-
logStatus(
468-
"[import_SILO] real zone range:"
469-
" [%d:%d, %d:%d, %d:%d] = %d x %d x %d zones",
470-
realZoneFirst.x,
471-
realZoneLast.x,
472-
realZoneFirst.y,
473-
realZoneLast.y,
474-
realZoneFirst.z,
475-
realZoneLast.z,
476-
dims.x,
477-
dims.y,
478-
dims.z);
479-
} else {
480-
dims = totalDims;
455+
if (hasGhostZones) {
456+
logStatus("[import_SILO] phony %s: [%d:%d, %d:%d, %d:%d]",
457+
isNodeCentered ? "nodes" : "cells",
458+
minIndex.x,
459+
maxIndex.x,
460+
minIndex.y,
461+
maxIndex.y,
462+
minIndex.z,
463+
maxIndex.z);
481464
}
482465

483-
if (dims.x < 1)
484-
dims.x = 1;
485-
if (dims.y < 1)
486-
dims.y = 1;
487-
if (dims.z < 1)
488-
dims.z = 1;
489-
490-
logStatus(
491-
"[import_SILO] total node dims:"
492-
" %d x %d x %d, total zone dims: %d x %d x %d, final dims: %d x %d x %d",
493-
totalNodeDims.x,
494-
totalNodeDims.y,
495-
totalNodeDims.z,
496-
totalDims.x,
497-
totalDims.y,
498-
totalDims.z,
499-
dims.x,
500-
dims.y,
501-
dims.z);
502-
503-
// Get origin and spacing from coordinates
466+
// Get origin and spacing from coordinates (use full extent including ghosts)
504467
float3 origin(0.f);
505468
float3 spacing(1.f);
469+
tsd::math::box3 roi;
506470

507-
// Spacing should be computed from consecutive real coordinates
508-
// Our grid might be rectilinear, let's compute some average spacing
509-
// that preserves the actual volume extent
471+
// Compute origin and spacing from the full grid extent
510472
for (int d = 0; d < 3; d++) {
511473
if (mesh->coords[d]) {
512474
if (mesh->datatype == DB_FLOAT) {
513475
float *c = (float *)mesh->coords[d];
514-
origin[d] = c[realNodeFirst[d]];
515-
if (mesh->dims[d] > realNodeFirst[d])
516-
spacing[d] = (c[realNodeLast[d]] - c[realNodeFirst[d]]) / dims[d];
476+
origin[d] = c[0]; // First node in full grid
477+
if (mesh->dims[d] > 1) {
478+
// Spacing from full grid extent
479+
spacing[d] = (c[mesh->dims[d] - 1] - c[0]) / nodes[d];
480+
}
481+
roi.lower[d] = c[minIndex[d]];
482+
roi.upper[d] = c[maxIndex[d]];
517483
} else if (mesh->datatype == DB_DOUBLE) {
518484
double *c = (double *)mesh->coords[d];
519-
origin[d] = c[realNodeFirst[d]];
520-
if (mesh->dims[d] > realNodeFirst[d])
521-
spacing[d] = (c[realNodeLast[d]] - c[realNodeFirst[d]]) / dims[d];
485+
origin[d] = c[0]; // First node in full grid
486+
if (mesh->dims[d] > 1) {
487+
// Spacing from full grid extent
488+
spacing[d] = (c[mesh->dims[d] - 1] - c[0]) / nodes[d];
489+
}
490+
roi.lower[d] = c[minIndex[d]];
491+
roi.upper[d] = c[maxIndex[d]];
522492
}
523493
}
524494
}
525495

526496
logStatus(
527-
"[import_SILO] extent: [%.9f, %.9f, %.9f ; %.9f, %.9f, %.9f] spacing: [%.9f, %.9f, %.9f]",
497+
"[import_SILO] origin: %.9f %.9f %.9f, spacing: %.9f %.9f %.9f, %s values",
528498
origin.x,
529499
origin.y,
530500
origin.z,
531-
origin.x + spacing.x * dims.x,
532-
origin.y + spacing.y * dims.y,
533-
origin.z + spacing.z * dims.z,
534501
spacing.x,
535502
spacing.y,
536-
spacing.z);
503+
spacing.z,
504+
isNodeCentered ? "nodes" : "cells");
537505

538506
field->setParameter("origin", origin);
539507
field->setParameter("spacing", spacing);
508+
field->setParameter("dataCentering", isNodeCentered ? "node" : "cell");
540509

541-
// Copy variable data - only copy non-ghost zones
542-
size_t numElements = dims.x * dims.y * dims.z;
543-
auto dataType = siloTypeToANARIType(var->datatype);
544-
auto dataArray = scene.createArray(dataType, dims.x, dims.y, dims.z);
545-
546-
void *dst = dataArray->map();
547-
510+
// Set ROI to exclude ghost zones if they exist
548511
if (hasGhostZones) {
549-
// Copy only the real zone sub-region defined by [realZoneFirst,
550-
// realZoneLast] For structured grids, min_index/max_index are NODE indices,
551-
// zones are nodes-1
552-
size_t dstIdx = 0;
553-
size_t elementSize = anari::sizeOf(dataType);
554-
555-
for (int k = realZoneFirst.z; k <= realZoneLast.z; k++) {
556-
for (int j = realZoneFirst.y; j <= realZoneLast.y; j++) {
557-
for (int i = realZoneFirst.x; i <= realZoneLast.x; i++) {
558-
int srcIdx = k * totalDims.y * totalDims.x + j * totalDims.x + i;
559-
std::memcpy((char *)dst + dstIdx * elementSize,
560-
(char *)var->vals[0] + srcIdx * elementSize,
561-
elementSize);
562-
dstIdx++;
563-
}
564-
}
565-
}
512+
field->setParameter("roi", roi);
566513
logStatus(
567-
"[import_SILO] copied %zu real zones (excluded ghosts from %d total)",
568-
dstIdx,
569-
totalDims.x * totalDims.y * totalDims.z);
570-
} else {
571-
// No ghost zones, direct copy
572-
void *src = var->vals[0];
573-
std::memcpy(dst, src, numElements * anari::sizeOf(dataType));
514+
"[import_SILO] ROI set to %.9f => %.9f; %.9f => %.9f; %.9f => %.9f to exclude ghost zones",
515+
roi.lower.x,
516+
roi.upper.x,
517+
roi.lower.y,
518+
roi.upper.y,
519+
roi.lower.z,
520+
roi.upper.z);
574521
}
575522

576-
dataArray->unmap();
577-
523+
// Load all variable data (including ghost zones)
524+
auto dataType = siloTypeToANARIType(var->datatype);
525+
ArrayRef dataArray;
526+
if (isNodeCentered) {
527+
dataArray = scene.createArray(dataType, nodes.x, nodes.y, nodes.z);
528+
} else {
529+
dataArray = scene.createArray(dataType, zones.x, zones.y, zones.z);
530+
}
531+
dataArray->setData(var->vals[0]);
578532
field->setParameterObject("data", *dataArray);
579533

580534
DBFreeQuadmesh(mesh);
@@ -876,6 +830,7 @@ static SpatialFieldRef import_SILO_multiMesh(Scene &scene,
876830
// Derived fields are always computed as float32
877831
anari::DataType globalDataType = ANARI_FLOAT32;
878832
bool firstBlock = true;
833+
bool isNodeCentered = true; // Default to node-centered
879834

880835
std::filesystem::path basePath = std::filesystem::path(file).parent_path();
881836

@@ -951,29 +906,35 @@ static SpatialFieldRef import_SILO_multiMesh(Scene &scene,
951906
if (qm->coords[0] && qm->dims[0] > realNodeFirst.x + 1) {
952907
if (qm->datatype == DB_FLOAT) {
953908
float *x = (float *)qm->coords[0];
954-
globalSpacing.x = x[realNodeFirst.x + 1] - x[realNodeFirst.x];
909+
globalSpacing.x = (x[realNodeLast.x] - x[realNodeFirst.x])
910+
/ (qm->dims[0] + (isNodeCentered ? -1 : 0));
955911
} else if (qm->datatype == DB_DOUBLE) {
956912
double *x = (double *)qm->coords[0];
957-
globalSpacing.x = x[realNodeFirst.x + 1] - x[realNodeFirst.x];
913+
globalSpacing.x = (x[realNodeLast.x] - x[realNodeFirst.x])
914+
/ (qm->dims[0] + (isNodeCentered ? -1 : 0));
958915
}
959916
}
960917
if (qm->coords[1] && qm->dims[1] > realNodeFirst.y + 1) {
961918
if (qm->datatype == DB_FLOAT) {
962919
float *y = (float *)qm->coords[1];
963-
globalSpacing.y = y[realNodeFirst.y + 1] - y[realNodeFirst.y];
920+
globalSpacing.y = (y[realNodeLast.y] - y[realNodeFirst.y])
921+
/ (qm->dims[1] + (isNodeCentered ? -1 : 0));
964922
} else if (qm->datatype == DB_DOUBLE) {
965923
double *y = (double *)qm->coords[1];
966-
globalSpacing.y = y[realNodeFirst.y + 1] - y[realNodeFirst.y];
924+
globalSpacing.y = (y[realNodeLast.y] - y[realNodeFirst.y])
925+
/ (qm->dims[1] + (isNodeCentered ? -1 : 0));
967926
}
968927
}
969928
if (qm->ndims > 2 && qm->coords[2]
970929
&& qm->dims[2] > realNodeFirst.z + 1) {
971930
if (qm->datatype == DB_FLOAT) {
972931
float *z = (float *)qm->coords[2];
973-
globalSpacing.z = z[realNodeFirst.z + 1] - z[realNodeFirst.z];
932+
globalSpacing.z = (z[realNodeLast.z] - z[realNodeFirst.z])
933+
/ (qm->dims[2] + (isNodeCentered ? -1 : 0));
974934
} else if (qm->datatype == DB_DOUBLE) {
975935
double *z = (double *)qm->coords[2];
976-
globalSpacing.z = z[realNodeFirst.z + 1] - z[realNodeFirst.z];
936+
globalSpacing.z = (z[realNodeLast.z] - z[realNodeFirst.z])
937+
/ (qm->dims[2] + (isNodeCentered ? -1 : 0));
977938
}
978939
}
979940
if (!isDerivedField) {
@@ -1007,6 +968,11 @@ static SpatialFieldRef import_SILO_multiMesh(Scene &scene,
1007968
info.origin = blockOrigin;
1008969
info.spacing = globalSpacing;
1009970

971+
// Capture centering from first block
972+
if (firstBlock) {
973+
isNodeCentered = (qv->centering == DB_NODECENT);
974+
}
975+
1010976
// Dims are in zones (nodes - 1)
1011977
info.dims = int3(realNodeLast.x - realNodeFirst.x,
1012978
realNodeLast.y - realNodeFirst.y,
@@ -1452,6 +1418,7 @@ static SpatialFieldRef import_SILO_multiMesh(Scene &scene,
14521418
field->setName(fieldName.c_str());
14531419
field->setParameter("origin", globalMin);
14541420
field->setParameter("spacing", globalSpacing);
1421+
field->setParameter("dataCentering", isNodeCentered ? "node" : "cell");
14551422
field->setParameterObject("data", *unifiedDataArray);
14561423

14571424
// Create single volume with unified field

0 commit comments

Comments
 (0)