vtk_io.h File Reference

#include "variables.h"
#include "logging.h"
#include "io.h"

Include dependency graph for vtk_io.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
PetscInt	CreateVTKFileFromMetadata (const char *filename, const VTKMetaData *meta, MPI_Comm comm)
	Creates and writes a VTK file (either .vts or .vtp) from a populated metadata struct.
PetscErrorCode	PrepareOutputCoordinates (UserCtx *user, PetscScalar **out_coords, PetscInt *out_nx, PetscInt *out_ny, PetscInt *out_nz, PetscInt *out_npoints)
	Creates a C array of coordinates corresponding to a subsampled (legacy-style) grid.
PetscErrorCode	PrepareOutputEulerianFieldData (UserCtx *user, Vec field_vec, PetscInt num_components, PetscScalar **out_data)
	Creates a C array of field data corresponding to a subsampled (legacy-style) grid.
PetscErrorCode	BeginStructuredVTKOutput (UserCtx *user, VTKMetaData *meta)
	Begins one structured VTK file: clears the metadata and builds its coordinates.
PetscErrorCode	AppendStructuredVTKField (UserCtx *user, VTKMetaData *meta, const char *name, Vec field_vec, PetscInt components)
	Adds one point-data field to a structured VTK file being assembled.
PetscErrorCode	FinishStructuredVTKOutput (VTKMetaData *meta, const char *filename)
	Writes an assembled structured VTK file and releases its buffers.
PetscErrorCode	PrepareOutputParticleData (UserCtx *user, PostProcessParams *pps, VTKMetaData *meta, PetscInt *p_n_total)
	Gathers, subsamples, and prepares all particle data for VTK output.

Function Documentation

CreateVTKFileFromMetadata()

PetscInt CreateVTKFileFromMetadata	(	const char *	filename,
		const VTKMetaData *	meta,
		MPI_Comm	comm
	)

Creates and writes a VTK file (either .vts or .vtp) from a populated metadata struct.

Parameters

[in]	filename	The output file name.
[in]	meta	Pointer to a VTKMetaData structure containing all necessary fields.
[in]	comm	The MPI communicator.

Returns

0 on success, non-zero on failure.

Creates and writes a VTK file (either .vts or .vtp) from a populated metadata struct.

Creates a VTK file from prepared metadata and field payloads.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/vtk_io.h.

See also

CreateVTKFileFromMetadata()

Definition at line 149 of file vtk_io.c.

{
    PetscMPIInt rank;
    MPI_Comm_rank(comm, &rank);
    PetscErrorCode ierr = 0;
 
    if (!rank) {
        char tmp_filename[PETSC_MAX_PATH_LEN];
        FILE *fp = NULL;
 
        ierr = PetscSNPrintf(tmp_filename, sizeof(tmp_filename), "%s.tmp", filename);CHKERRQ(ierr);
        LOG_ALLOW(GLOBAL, LOG_INFO, "Rank 0 writing combined VTK file '%s'.\n", filename);
        fp = fopen(tmp_filename, "wb");
        if (!fp) {
            LOG_ALLOW(GLOBAL, LOG_ERROR, "fopen failed for %s.\n", tmp_filename);
            return PETSC_ERR_FILE_OPEN;
        }
 
        PetscInt boffset = 0;
 
        ierr = WriteVTKFileHeader(fp, meta, &boffset);
        if (ierr) {
            fclose(fp);
            remove(tmp_filename);
            return ierr;
        }
 
        if (meta->coords) {
            ierr = WriteVTKAppendedBlock(fp, meta->coords, 3 * meta->npoints, sizeof(PetscScalar));
            if (ierr) {
                fclose(fp);
                remove(tmp_filename);
                return ierr;
            }
        }
 
        for (PetscInt i = 0; i < meta->num_point_data_fields; i++) {
            const VTKFieldInfo *field = &meta->point_data_fields[i];
            if (field->data) {
                ierr = WriteVTKAppendedBlock(fp, field->data, field->num_components * meta->npoints, sizeof(PetscScalar));
                if (ierr) {
                    fclose(fp);
                    remove(tmp_filename);
                    return ierr;
                }
            }
        }
        if (meta->fileType == VTK_POLYDATA) {
            if (meta->connectivity) {
                ierr = WriteVTKAppendedBlock(fp, meta->connectivity, meta->npoints, sizeof(PetscInt));
                if (ierr) {
                    fclose(fp);
                    remove(tmp_filename);
                    return ierr;
                }
            }
            if (meta->offsets) {
                ierr = WriteVTKAppendedBlock(fp, meta->offsets, meta->npoints, sizeof(PetscInt));
                if (ierr) {
                    fclose(fp);
                    remove(tmp_filename);
                    return ierr;
                }
            }
        }
        ierr = WriteVTKFileFooter(fp, meta);
        if (ierr) {
            fclose(fp);
            remove(tmp_filename);
            return ierr;
        }
        if (fclose(fp) != 0) {
            remove(tmp_filename);
            return PETSC_ERR_FILE_WRITE;
        }
        if (rename(tmp_filename, filename) != 0) {
            remove(tmp_filename);
            return PETSC_ERR_FILE_WRITE;
        }
        LOG_ALLOW(GLOBAL, LOG_INFO, "Rank 0 finished writing VTK file '%s'.\n", filename);
    }
    return 0;
}

Here is the caller graph for this function:

PrepareOutputCoordinates()

PetscErrorCode PrepareOutputCoordinates	(	UserCtx *	user,
		PetscScalar **	out_coords,
		PetscInt *	out_nx,
		PetscInt *	out_ny,
		PetscInt *	out_nz,
		PetscInt *	out_npoints
	)

Creates a C array of coordinates corresponding to a subsampled (legacy-style) grid.

This function gathers the full, distributed grid coordinates onto rank 0. On rank 0, it then allocates a new, smaller C array and copies only the coordinates for the nodes within the range [0..IM-2, 0..JM-2, 0..KM-2]. This produces a contiguous array of points for a grid of size (IM-1)x(JM-1)x(KM-1), matching the legacy output. The output arrays are only allocated and valid on rank 0.

Parameters

[in]	user	The UserCtx containing the grid information (DM, IM/JM/KM).
[out]	out_coords	On rank 0, a pointer to the newly allocated C array for coordinate data. NULL on other ranks.
[out]	out_nx	The number of points in the x-dimension for the new grid (IM-1).
[out]	out_ny	The number of points in the y-dimension for the new grid (JM-1).
[out]	out_nz	The number of points in the z-dimension for the new grid (KM-1).
[out]	out_npoints	The total number of points in the new grid.

Returns

PetscErrorCode

Creates a C array of coordinates corresponding to a subsampled (legacy-style) grid.

Local to this translation unit.

Definition at line 241 of file vtk_io.c.

{
    PetscErrorCode ierr;
    Vec            coords_global;
    PetscInt       IM, JM, KM;
    PetscInt       N_full_coords;
    PetscScalar    *full_coords_arr = NULL;
 
    PetscFunctionBeginUser;
    ierr = DMDAGetInfo(user->da, NULL, &IM, &JM, &KM, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); CHKERRQ(ierr);
 
    // Set the dimensions of the output grid
    *out_nx = IM - 1;
    *out_ny = JM - 1;
    *out_nz = KM - 1;
    *out_npoints = (*out_nx) * (*out_ny) * (*out_nz);
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Preparing subsampled coordinates for a %" PetscInt_FMT "x%" PetscInt_FMT "x%" PetscInt_FMT " grid.\n", *out_nx, *out_ny, *out_nz);
 
    // --- Step 1: Gather the full coordinate vector to rank 0 ---
    ierr = DMGetCoordinates(user->da, &coords_global); CHKERRQ(ierr);
    // Reuse of your existing, proven utility function's logic
    ierr = VecToArrayOnRank0(coords_global, &N_full_coords, &full_coords_arr); CHKERRQ(ierr);
 
    // --- Step 2: On rank 0, subsample the gathered array ---
    if (user->simCtx->rank == 0) {
        // We expect N_full_coords to be IM * JM * KM * 3
        if (N_full_coords != IM * JM * KM * 3) {
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED, "Gathered coordinate array has wrong size. Expected %" PetscInt_FMT ", got %" PetscInt_FMT, IM * JM * KM * 3, N_full_coords);
        }
 
        // Allocate the smaller output C array
        ierr = PetscMalloc1(3 * (*out_npoints), out_coords); CHKERRQ(ierr);
 
        // Loop over the smaller grid dimensions and copy data
        PetscInt p_out = 0; // Index for the small output array
        for (PetscInt k = 0; k < *out_nz; k++) {
            for (PetscInt j = 0; j < *out_ny; j++) {
                for (PetscInt i = 0; i < *out_nx; i++) {
                    // Calculate the index in the full, 1D source array
                    PetscInt p_in = 3 * (k * (JM * IM) + j * IM + i);
                    
                    (*out_coords)[p_out++] = full_coords_arr[p_in + 0]; // x
                    (*out_coords)[p_out++] = full_coords_arr[p_in + 1]; // y
                    (*out_coords)[p_out++] = full_coords_arr[p_in + 2]; // z
                }
            }
        }
        // Free the temporary full array that was allocated by VecToArrayOnRank0
        ierr = PetscFree(full_coords_arr); CHKERRQ(ierr);
    } else {
        // On other ranks, the output pointer is NULL
        *out_coords = NULL;
    }
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Subsampled coordinates prepared on rank 0 with total %" PetscInt_FMT " points.\n", *out_npoints);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PrepareOutputEulerianFieldData()

PetscErrorCode PrepareOutputEulerianFieldData	(	UserCtx *	user,
		Vec	field_vec,
		PetscInt	num_components,
		PetscScalar **	out_data
	)

Creates a C array of field data corresponding to a subsampled (legacy-style) grid.

This function gathers a full, distributed PETSc vector to rank 0. On rank 0, it then allocates a new, smaller C array and copies only the data components for nodes within the range [0..IM-2, 0..JM-2, 0..KM-2]. This produces a contiguous data array that perfectly matches the point ordering of the subsampled coordinates. The output array is only allocated and valid on rank 0.

Parameters

[in]	user	The UserCtx for grid information.
[in]	field_vec	The full-sized PETSc vector containing the field data (e.g., user->P_nodal).
[in]	num_components	The number of components for this field (1 for scalar, 3 for vector).
[out]	out_data	On rank 0, a pointer to the newly allocated C array for the field data. NULL on other ranks.

Returns

PetscErrorCode

Creates a C array of field data corresponding to a subsampled (legacy-style) grid.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/vtk_io.h.

See also

PrepareOutputEulerianFieldData()

Definition at line 377 of file vtk_io.c.

{
    PetscErrorCode ierr;
    PetscInt       IM, JM, KM, nx, ny, nz, npoints;
    PetscInt       N_full_field = 0;
    PetscScalar   *full_field_arr = NULL;
 
    PetscFunctionBeginUser;
    ierr = DMDAGetInfo(user->da, NULL, &IM, &JM, &KM, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL); CHKERRQ(ierr);
    nx = IM - 1; ny = JM - 1; nz = KM - 1;
    npoints = nx * ny * nz;
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG,
              "Preparing subsampled field data (dof=%" PetscInt_FMT ") for a %" PetscInt_FMT "x%" PetscInt_FMT "x%" PetscInt_FMT " grid.\n",
              num_components, nx, ny, nz);
 
    // --- Step 1: Gather the full field vector to rank 0 ---
    ierr = VecToArrayOnRank0(field_vec, &N_full_field, &full_field_arr); CHKERRQ(ierr);
 
    // --- Step 2: On rank 0, subsample the gathered array ---
    if (user->simCtx->rank == 0) {
        // Sanity check the size of the gathered array
        if (N_full_field != IM * JM * KM * num_components) {
            SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_UNEXPECTED,
                    "Gathered field array has wrong size. Expected %" PetscInt_FMT ", got %" PetscInt_FMT,
                    IM * JM * KM * num_components, N_full_field);
        }
 
        // Allocate the smaller output C array (AoS layout expected by your writer)
        ierr = PetscMalloc1(num_components * npoints, out_data); CHKERRQ(ierr);
 
        /*
        // ======== Layout diagnostics (rank 0 only, read-only) ========
        if (num_components > 1 && full_field_arr) {
            const PetscInt n_full = IM * JM * KM;
 
            // Choose a safe interior probe point and its +i neighbor
            const PetscInt ic = (IM > 4) ? IM/2 : 1;
            const PetscInt jc = (JM > 4) ? JM/2 : 1;
            const PetscInt kc = (KM > 4) ? KM/2 : 1;
            const PetscInt ic1 = (ic + 1 < IM-1) ? (ic + 1) : ic; // stay interior
 
            const PetscInt t_full_c  = (kc * JM + jc) * IM + ic;
            const PetscInt t_full_c1 = (kc * JM + jc) * IM + ic1;
 
            // If source were AoS:  full[num_components * t + comp]
            PetscScalar aos_c0   = full_field_arr[num_components * t_full_c  + 0];
            PetscScalar aos_c1   = full_field_arr[num_components * t_full_c1 + 0];
            PetscScalar aos_c0_y = (num_components > 1) ? full_field_arr[num_components * t_full_c + 1] : 0.0;
            PetscScalar aos_c0_z = (num_components > 2) ? full_field_arr[num_components * t_full_c + 2] : 0.0;
 
            // If source were SoA:  full[comp * n_full + t]
            PetscScalar soa_c0   = full_field_arr[0 * n_full + t_full_c];
            PetscScalar soa_c1   = full_field_arr[0 * n_full + t_full_c1];
            PetscScalar soa_c0_y = (num_components > 1) ? full_field_arr[1 * n_full + t_full_c] : 0.0;
            PetscScalar soa_c0_z = (num_components > 2) ? full_field_arr[2 * n_full + t_full_c] : 0.0;
 
            double dAOS = fabs((double)(aos_c1 - aos_c0)); // spatial delta along +i
            double dSOA = fabs((double)(soa_c1 - soa_c0));
 
            const char *verdict = (dSOA < dAOS) ? "LIKELY SoA (component-major)" : "LIKELY AoS (interleaved)";
 
            LOG_ALLOW(GLOBAL, LOG_INFO,
                "Layout probe @ (i=%" PetscInt_FMT ", j=%" PetscInt_FMT ", k=%" PetscInt_FMT ") vs +i neighbor:\n"
                "  AoS-candidate:  Ux(center)=%.6e  Ux(+i)=%.6e  |Δ|=%.6e ; Uy(center)=%.6e  Uz(center)=%.6e\n"
                "  SoA-candidate:  Ux(center)=%.6e  Ux(+i)=%.6e  |Δ|=%.6e ; Uy(center)=%.6e  Uz(center)=%.6e\n"
                "  Verdict: %s\n",
                ic, jc, kc,
                (double)aos_c0,  (double)aos_c1,  dAOS, (double)aos_c0_y, (double)aos_c0_z,
                (double)soa_c0,  (double)soa_c1,  dSOA, (double)soa_c0_y, (double)soa_c0_z,
                verdict
            );
        } else if (num_components == 1) {
            LOG_ALLOW(GLOBAL, LOG_INFO, "Layout probe: scalar field (dof=1) — no component interleave to detect.\n");
        }
        // ======== END diagnostics ========
        */
 
        // --- PACK: assumes source is already AoS ---
        {
            PetscInt p_out = 0;
            for (PetscInt k = 0; k < nz; k++) {
                for (PetscInt j = 0; j < ny; j++) {
                    for (PetscInt i = 0; i < nx; i++) {
                        const PetscInt p_in_start = num_components * (k * (JM * IM) + j * IM + i);
                        for (PetscInt c = 0; c < num_components; c++) {
                            (*out_data)[p_out++] = full_field_arr[p_in_start + c];
                        }
                    }
                }
            }
        }
 
        // Free temporary gathered array
        ierr = PetscFree(full_field_arr); CHKERRQ(ierr);
    } else {
        // Other ranks: no allocation
        *out_data = NULL;
    }
 
    LOG_ALLOW(GLOBAL, LOG_DEBUG, "Subsampled field data prepared on rank 0 with total %" PetscInt_FMT " points.\n", npoints);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

BeginStructuredVTKOutput()

PetscErrorCode BeginStructuredVTKOutput	(	UserCtx *	user,
		VTKMetaData *	meta
	)

Begins one structured VTK file: clears the metadata and builds its coordinates.

Every Eulerian VTK file is assembled the same way — coordinates, then a set of point-data fields, then the write. This trio exists so the two producers share that shape instead of restating it, and so the buffers each producer allocates are released in one place rather than in neither.

Parameters

[in]	user	Block context supplying the grid.
[out]	meta	Metadata to initialize.

Returns

Zero on success, or a PETSc error.

Begins one structured VTK file: clears the metadata and builds its coordinates.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/vtk_io.h.

See also

BeginStructuredVTKOutput()

Definition at line 308 of file vtk_io.c.

{
    PetscFunctionBeginUser;
    PetscCheck(user != NULL && meta != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Context and metadata are required.");
    PetscCall(PetscMemzero(meta, sizeof(*meta)));
    meta->fileType = VTK_STRUCTURED;
    PetscCall(PrepareOutputCoordinates(user, &meta->coords, &meta->mx, &meta->my, &meta->mz,
                                       &meta->npoints));
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

AppendStructuredVTKField()

PetscErrorCode AppendStructuredVTKField	(	UserCtx *	user,
		VTKMetaData *	meta,
		const char *	name,
		Vec	field_vec,
		PetscInt	components
	)

Adds one point-data field to a structured VTK file being assembled.

Skips the field with a warning when the file is already at its field limit, so a long output list degrades rather than corrupting the metadata.

Parameters

[in]	user	Block context supplying the grid.
[in,out]	meta	Metadata to append to.
[in]	name	Field name as it will appear in the file.
[in]	field_vec	Nodal vector holding the values.
[in]	components	One for a scalar, three for a vector.

Returns

Zero on success, or a PETSc error.

Adds one point-data field to a structured VTK file being assembled.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/vtk_io.h.

See also

AppendStructuredVTKField()

Definition at line 326 of file vtk_io.c.

{
    VTKFieldInfo *entry = NULL;
 
    PetscFunctionBeginUser;
    PetscCheck(user != NULL && meta != NULL && name != NULL && field_vec != NULL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "Context, metadata, name, and vector are required.");
    if (meta->num_point_data_fields >= MAX_POINT_DATA_FIELDS) {
        LOG_ALLOW(GLOBAL, LOG_WARNING,
                  "Reached the %d point-data field limit; '%s' is not written.\n",
                  MAX_POINT_DATA_FIELDS, name);
        PetscFunctionReturn(0);
    }
    entry = &meta->point_data_fields[meta->num_point_data_fields];
    PetscCall(PetscStrncpy(entry->name, name, MAX_VTK_FIELD_NAME_LENGTH));
    entry->num_components = components;
    PetscCall(PrepareOutputEulerianFieldData(user, field_vec, components, &entry->data));
    meta->num_point_data_fields++;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

FinishStructuredVTKOutput()

PetscErrorCode FinishStructuredVTKOutput	(	VTKMetaData *	meta,
		const char *	filename
	)

Writes an assembled structured VTK file and releases its buffers.

Parameters

[in,out]	meta	Metadata to write and then release.
[in]	filename	Destination path.

Returns

Zero on success, or a PETSc error.

Writes an assembled structured VTK file and releases its buffers.

Full API contract (arguments, ownership, side effects) is documented with the header declaration in include/vtk_io.h.

See also

FinishStructuredVTKOutput()

Definition at line 354 of file vtk_io.c.

{
    PetscFunctionBeginUser;
    PetscCheck(meta != NULL && filename != NULL, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL,
               "Metadata and filename are required.");
    PetscCall(CreateVTKFileFromMetadata(filename, meta, PETSC_COMM_WORLD));
    for (PetscInt index = 0; index < meta->num_point_data_fields; ++index) {
        if (meta->point_data_fields[index].data) {
            PetscCall(PetscFree(meta->point_data_fields[index].data));
        }
    }
    if (meta->coords) PetscCall(PetscFree(meta->coords));
    meta->num_point_data_fields = 0;
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function:

PrepareOutputParticleData()

PetscErrorCode PrepareOutputParticleData	(	UserCtx *	user,
		PostProcessParams *	pps,
		VTKMetaData *	meta,
		PetscInt *	p_n_total
	)

Gathers, subsamples, and prepares all particle data for VTK output.

This function is a COLLECTIVE operation. All ranks must enter it. The heavy lifting (memory allocation, subsampling) is performed only on rank 0.

1. Gathers the full coordinate and field data from the distributed DMSwarm to rank 0.
2. On rank 0, subsamples the data based on pps->particle_output_freq.
3. On rank 0, populates the VTKMetaData struct with the new, smaller, subsampled data arrays.

Parameters

[in]	user	The UserCtx containing the DMSwarm.
[in]	pps	The PostProcessParams struct for configuration.
[out]	meta	A pointer to the VTKMetaData struct to be populated (on rank 0).
[out]	p_n_total	On rank 0, the total number of particles before subsampling.

Returns

PetscErrorCode

Gathers, subsamples, and prepares all particle data for VTK output.

Local to this translation unit.

Definition at line 486 of file vtk_io.c.

{
    PetscErrorCode ierr;
    PetscInt n_total_particles, n_components;
    PetscInt stride = pps->particle_output_freq > 0 ? pps->particle_output_freq : 1;
 
    PetscFunctionBeginUser;
 
    // Initialize output parameters on all ranks
    *p_n_total = 0;
 
    // --- The entire preparation process is a series of collective gathers followed by rank-0 processing ---
 
    // --- Step 1: Gather Full Coordinates from SOURCE Swarm (Collective) ---
    // This establishes the ground truth for particle positions and total count.
    PetscReal *full_coords_arr = NULL;
    PetscDataType coordinate_type;
    // NOTE: This assumes SwarmFieldToArrayOnRank0 exists and works for PetscScalar fields.
    // If not, it may need to be generalized like the logic below.
    ierr = SwarmFieldToArrayOnRank0(user->swarm, ParticleFieldName(PARTICLE_FIELD_ID_POSITION),
                                    &n_total_particles, &n_components, &coordinate_type,
                                    (void**)&full_coords_arr); CHKERRQ(ierr);
    PetscCheck(coordinate_type == PETSC_REAL,
               PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
               "Particle position field must use PETSC_REAL storage, not %s.",
               PetscDataTypes[coordinate_type]);
 
    *p_n_total = n_total_particles;
    if (n_total_particles == 0) {
        ierr = PetscFree(full_coords_arr); CHKERRQ(ierr);
        PetscFunctionReturn(0);
    }
    PetscCheck(n_components == 3,
               PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG,
               "Coordinate field position must have 3 components, but has %" PetscInt_FMT,
               n_components);
 
    // --- Step 2: Prepare and Subsample Coordinates (Rank 0 Only) ---
    if (user->simCtx->rank == 0) {
        // --- Subsampling Calculation ---
        meta->npoints = (n_total_particles > 0) ? (n_total_particles - 1) / stride + 1 : 0;
 
        LOG_ALLOW(LOCAL, LOG_DEBUG, "Subsampling %" PetscInt_FMT " total particles with stride %" PetscInt_FMT " -> %" PetscInt_FMT " output particles.\n",
                  n_total_particles, stride, meta->npoints);
 
        // --- Prepare Final Coordinates Array ---
        ierr = PetscMalloc1(3 * meta->npoints, &meta->coords); CHKERRQ(ierr);
        for (PetscInt i = 0; i < meta->npoints; i++) {
            PetscInt source_idx = i * stride;
            for (int d = 0; d < 3; d++) meta->coords[3 * i + d] = full_coords_arr[3 * source_idx + d];
        }
    }
    ierr = PetscFree(full_coords_arr); CHKERRQ(ierr);
 
    // --- Step 3: Gather Fields (Collective), then Subsample and Store (Rank 0 Only) ---
    char *fields_copy, *field_name;
    PetscInt num_fields = 0;
    ierr = PetscStrallocpy(pps->particle_fields, &fields_copy); CHKERRQ(ierr);
    field_name = strtok(fields_copy, ",");
    while (field_name && num_fields < MAX_POINT_DATA_FIELDS) {
        TrimWhitespace(field_name);
        if (!*field_name || strcasecmp(field_name, "position") == 0) {
            field_name = strtok(NULL, ","); continue;
        }
 
        // A. Determine which swarm is the source for this field on every rank.
        DM swarm_to_use = user->swarm; // Default to the main solver swarm.
        const char* internal_name = field_name; // Default internal name to user-facing name.
 
        PetscErrorCode check_ierr;
        ierr = PetscPushErrorHandler(PetscIgnoreErrorHandler, NULL); CHKERRQ(ierr);
        check_ierr = DMSwarmGetField(user->post_swarm, field_name, NULL, NULL, NULL);
        ierr = PetscPopErrorHandler(); CHKERRQ(ierr);
        if (!check_ierr) {
            swarm_to_use = user->post_swarm;
            ierr = DMSwarmRestoreField(user->post_swarm, field_name, NULL, NULL, NULL); CHKERRQ(ierr);
            if (user->simCtx->rank == 0) {
                LOG_ALLOW(LOCAL, LOG_DEBUG, "Field '%s' will be sourced from the post-processing swarm.\n", field_name);
            }
        } else {
            if (user->simCtx->rank == 0) {
                LOG_ALLOW(LOCAL, LOG_DEBUG, "Field '%s' will be sourced from the main solver swarm.\n", field_name);
            }
            if (strcasecmp(field_name, "pid") == 0) internal_name = "DMSwarm_pid";
            else if (strcasecmp(field_name, "CellID") == 0) internal_name = "DMSwarm_CellID";
            else if (strcasecmp(field_name, "Migration Status") == 0) internal_name = "DMSwarm_location_status";
        }
 
        // B. Gather the complete field on rank 0 before applying the output stride.
        // SwarmFieldToArrayOnRank0 handles both scalar and vector field layouts.
        void* full_field_arr_void = NULL;
        PetscInt field_total_particles, field_num_components;
        PetscDataType field_data_type;
 
        ierr = SwarmFieldToArrayOnRank0(swarm_to_use, internal_name,
                                        &field_total_particles, &field_num_components,
                                        &field_data_type, &full_field_arr_void); CHKERRQ(ierr);
 
        if (field_total_particles != n_total_particles) {
            if (user->simCtx->rank == 0) {
                LOG_ALLOW(LOCAL, LOG_WARNING, "Field '%s' has %" PetscInt_FMT " particles, but expected %" PetscInt_FMT ". Skipping.\n", field_name, field_total_particles, n_total_particles);
            }
            ierr = PetscFree(full_field_arr_void); CHKERRQ(ierr);
            field_name = strtok(NULL, ","); continue;
        }
 
        if (user->simCtx->rank == 0) {
            // C. Allocate final array (ALWAYS as PetscScalar) and copy/cast subsampled data.
            VTKFieldInfo* current_field = &meta->point_data_fields[meta->num_point_data_fields];
            strncpy(current_field->name, field_name, MAX_VTK_FIELD_NAME_LENGTH - 1);
            current_field->name[MAX_VTK_FIELD_NAME_LENGTH - 1] = '\0';
            current_field->num_components = field_num_components;
            ierr = PetscMalloc1(current_field->num_components * meta->npoints, (PetscScalar**)&current_field->data); CHKERRQ(ierr);
 
            PetscScalar* final_data_arr = (PetscScalar*)current_field->data;
 
            // D. Perform the subsampling and potential casting
            if (field_data_type == PETSC_INT64) {
                PetscInt64* source_arr = (PetscInt64*)full_field_arr_void;
                for (PetscInt i = 0; i < meta->npoints; i++) {
                    final_data_arr[i] = (PetscScalar)source_arr[i * stride];
                }
            } else if (field_data_type == PETSC_INT) {
                PetscInt* source_arr = (PetscInt*)full_field_arr_void;
                for (PetscInt i = 0; i < meta->npoints; i++) {
                    PetscInt source_idx = i * stride;
                    for (PetscInt c = 0; c < current_field->num_components; c++) {
                        final_data_arr[current_field->num_components * i + c] = (PetscScalar)source_arr[current_field->num_components * source_idx + c];
                    }
                }
            } else if (field_data_type == PETSC_REAL) {
                PetscReal* source_arr = (PetscReal*)full_field_arr_void;
                for (PetscInt i = 0; i < meta->npoints; i++) {
                    PetscInt source_idx = i * stride;
                    for (PetscInt c = 0; c < current_field->num_components; c++) {
                        final_data_arr[current_field->num_components * i + c] = source_arr[current_field->num_components * source_idx + c];
                    }
                }
#if defined(PETSC_USE_COMPLEX)
            } else if (field_data_type == PETSC_SCALAR) {
                PetscScalar* source_arr = (PetscScalar*)full_field_arr_void;
                for (PetscInt i = 0; i < meta->npoints; i++) {
                    PetscInt source_idx = i * stride;
                    for (PetscInt c = 0; c < current_field->num_components; c++) {
                        final_data_arr[current_field->num_components * i + c] = source_arr[current_field->num_components * source_idx + c];
                    }
                }
#endif
            } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP,
                           "VTK particle output does not support field '%s' with PETSc type %s.",
                           field_name, PetscDataTypes[field_data_type]);
 
            meta->num_point_data_fields++;
        }
        ierr = PetscFree(full_field_arr_void); CHKERRQ(ierr);
        num_fields++;
        field_name = strtok(NULL, ",");
    }
    ierr = PetscFree(fields_copy); CHKERRQ(ierr);
 
    // --- Step 4: Finalize VTK MetaData for PolyData (Rank 0 Only) ---
    if (user->simCtx->rank == 0) {
        if (meta->npoints > 0) {
            meta->fileType = VTK_POLYDATA;
            ierr = PetscMalloc1(meta->npoints, &meta->connectivity); CHKERRQ(ierr);
            ierr = PetscMalloc1(meta->npoints, &meta->offsets); CHKERRQ(ierr);
            for (PetscInt i = 0; i < meta->npoints; i++) {
                meta->connectivity[i] = i;
                meta->offsets[i] = i + 1;
            }
        }
    }
 
    ierr = MPI_Barrier(PETSC_COMM_WORLD); CHKERRQ(ierr);
 
    PetscFunctionReturn(0);
}

Here is the call graph for this function:

Here is the caller graph for this function: