Line data Source code
1 : #include "cfd/io/vtk_output.h"
2 : #include "cfd/core/filesystem.h"
3 : #include "cfd/core/grid.h"
4 : #include "cfd/core/indexing.h"
5 : #include "cfd/core/logging.h"
6 :
7 :
8 : #include "cfd/solvers/navier_stokes_solver.h"
9 : #include "vtk_output_internal.h"
10 : #include <stdio.h>
11 : #include <stdlib.h>
12 : #include <string.h>
13 :
14 :
15 : //=============================================================================
16 : // UTILITY FUNCTIONS
17 : //=============================================================================
18 :
19 : // Helper to build filepath
20 2 : static void build_filepath(char* filepath, size_t filepath_size, const char* run_dir,
21 : const char* filename) {
22 : #ifdef _WIN32
23 : snprintf(filepath, filepath_size, "%s\\%s", run_dir, filename);
24 : #else
25 2 : snprintf(filepath, filepath_size, "%s/%s", run_dir, filename);
26 : #endif
27 2 : }
28 :
29 : //=============================================================================
30 : // VTK OUTPUT DISPATCH
31 : //=============================================================================
32 :
33 : // Function pointer type for VTK output handlers
34 : typedef void (*vtk_output_handler)(const char* run_dir, const char* prefix, int step,
35 : const flow_field* field, const grid* grid);
36 :
37 : // VTK output handler functions
38 1 : static void write_velocity_vtk(const char* run_dir, const char* prefix, int step,
39 : const flow_field* field, const grid* grid) {
40 1 : const char* name = prefix ? prefix : "velocity";
41 1 : char filename[256], filepath[512];
42 :
43 1 : snprintf(filename, sizeof(filename), "%s_%03d.vtk", name, step);
44 1 : build_filepath(filepath, sizeof(filepath), run_dir, filename);
45 :
46 1 : write_vtk_vector_output(filepath, "velocity", field->u, field->v, field->w, grid->nx, grid->ny,
47 1 : grid->nz, grid->xmin, grid->xmax, grid->ymin, grid->ymax, grid->zmin,
48 1 : grid->zmax);
49 1 : }
50 :
51 0 : static void write_full_field_vtk(const char* run_dir, const char* prefix, int step,
52 : const flow_field* field, const grid* grid) {
53 0 : const char* name = prefix ? prefix : "flow_field";
54 0 : char filename[256], filepath[512];
55 :
56 0 : snprintf(filename, sizeof(filename), "%s_%03d.vtk", name, step);
57 0 : build_filepath(filepath, sizeof(filepath), run_dir, filename);
58 :
59 0 : write_vtk_flow_field(filepath, field, grid->nx, grid->ny, grid->nz, grid->xmin, grid->xmax,
60 0 : grid->ymin, grid->ymax, grid->zmin, grid->zmax);
61 0 : }
62 :
63 : // VTK output handler table - indexed by VtkOutputType
64 : // Note: VTK_OUTPUT_VELOCITY_MAGNITUDE is handled via vtk_write_scalar_field, not dispatch
65 : static const vtk_output_handler vtk_output_table[] = {
66 : NULL, // VTK_OUTPUT_VELOCITY_MAGNITUDE = 0 (uses vtk_write_scalar_field)
67 : write_velocity_vtk, // VTK_OUTPUT_VELOCITY = 1
68 : write_full_field_vtk // VTK_OUTPUT_FULL_FIELD = 2
69 : };
70 :
71 : #define VTK_OUTPUT_TABLE_SIZE (sizeof(vtk_output_table) / sizeof(vtk_output_table[0]))
72 :
73 1 : void vtk_dispatch_output(vtk_output_type vtk_type, const char* run_dir, const char* prefix,
74 : int step, const flow_field* field, const grid* grid) {
75 : // VTK_OUTPUT_VELOCITY_MAGNITUDE should use vtk_write_scalar_field directly
76 1 : if (vtk_type == VTK_OUTPUT_VELOCITY_MAGNITUDE) {
77 0 : cfd_warning("VTK_OUTPUT_VELOCITY_MAGNITUDE should use vtk_write_scalar_field");
78 0 : return;
79 : }
80 :
81 : // Bounds check and dispatch via function table
82 1 : if (vtk_type < VTK_OUTPUT_TABLE_SIZE && vtk_output_table[vtk_type] != NULL) {
83 1 : vtk_output_table[vtk_type](run_dir, prefix, step, field, grid);
84 : } else {
85 0 : cfd_warning("Unknown VTK output type");
86 : }
87 : }
88 :
89 : // Write pre-computed scalar field to VTK
90 1 : void vtk_write_scalar_field(const char* run_dir, const char* prefix, int step,
91 : const char* field_name, const double* data, const grid* grid) {
92 1 : if (!run_dir || !data || !grid) {
93 0 : return;
94 : }
95 :
96 1 : const char* name = prefix ? prefix : "scalar";
97 1 : char filename[256], filepath[512];
98 :
99 1 : snprintf(filename, sizeof(filename), "%s_%03d.vtk", name, step);
100 1 : build_filepath(filepath, sizeof(filepath), run_dir, filename);
101 :
102 1 : write_vtk_output(filepath, field_name, data, grid->nx, grid->ny, grid->nz, grid->xmin,
103 1 : grid->xmax, grid->ymin, grid->ymax, grid->zmin, grid->zmax);
104 : }
105 :
106 : //=============================================================================
107 : // VTK LEGACY API
108 : //=============================================================================
109 :
110 12 : void write_vtk_output(const char* filename, const char* field_name, const double* data, size_t nx,
111 : size_t ny, size_t nz, double xmin, double xmax, double ymin, double ymax,
112 : double zmin, double zmax) {
113 12 : if (!filename || !field_name || !data || nx < 2 || ny < 2 || nz < 1 ||
114 9 : xmax <= xmin || ymax <= ymin || (nz > 1 && zmax <= zmin)) {
115 : return;
116 : }
117 :
118 9 : FILE* fp = fopen(filename, "w");
119 9 : if (fp == NULL) {
120 0 : cfd_error("Failed to open VTK output file");
121 0 : return;
122 : }
123 :
124 9 : double dz = (nz > 1) ? (zmax - zmin) / (double)(nz - 1) : 1.0;
125 :
126 : // Write VTK header
127 9 : fprintf(fp, "# vtk DataFile Version 3.0\n");
128 9 : fprintf(fp, "CFD Framework Output\n");
129 9 : fprintf(fp, "ASCII\n");
130 9 : fprintf(fp, "DATASET STRUCTURED_POINTS\n");
131 9 : fprintf(fp, "DIMENSIONS %zu %zu %zu\n", nx, ny, nz);
132 9 : fprintf(fp, "ORIGIN %f %f %f\n", xmin, ymin, zmin);
133 9 : fprintf(fp, "SPACING %f %f %f\n", (xmax - xmin) / (nx - 1), (ymax - ymin) / (ny - 1), dz);
134 :
135 :
136 : // Write field data
137 9 : fprintf(fp, "\nPOINT_DATA %zu\n", nx * ny * nz);
138 9 : fprintf(fp, "SCALARS %s float 1\n", field_name);
139 9 : fprintf(fp, "LOOKUP_TABLE default\n");
140 :
141 19 : for (size_t k = 0; k < nz; k++) {
142 74 : for (size_t j = 0; j < ny; j++) {
143 584 : for (size_t i = 0; i < nx; i++) {
144 520 : size_t idx = (k * nx * ny) + IDX_2D(i, j, nx);
145 520 : fprintf(fp, "%f\n", data[idx]);
146 : }
147 : }
148 : }
149 :
150 9 : fclose(fp);
151 : }
152 :
153 9 : void write_vtk_vector_output(const char* filename, const char* field_name, const double* u_data,
154 : const double* v_data, const double* w_data, size_t nx, size_t ny,
155 : size_t nz, double xmin, double xmax, double ymin, double ymax,
156 : double zmin, double zmax) {
157 9 : if (!filename || !field_name || !u_data || !v_data || nx < 2 || ny < 2 || nz < 1 ||
158 5 : xmax <= xmin || ymax <= ymin || (nz > 1 && zmax <= zmin)) {
159 : return;
160 : }
161 :
162 5 : FILE* fp = fopen(filename, "w");
163 5 : if (fp == NULL) {
164 0 : cfd_error("Failed to open VTK vector output file");
165 0 : return;
166 : }
167 :
168 5 : double dz = (nz > 1) ? (zmax - zmin) / (double)(nz - 1) : 1.0;
169 :
170 : // Write VTK header
171 5 : fprintf(fp, "# vtk DataFile Version 3.0\n");
172 5 : fprintf(fp, "CFD Framework Vector Output\n");
173 5 : fprintf(fp, "ASCII\n");
174 5 : fprintf(fp, "DATASET STRUCTURED_POINTS\n");
175 5 : fprintf(fp, "DIMENSIONS %zu %zu %zu\n", nx, ny, nz);
176 5 : fprintf(fp, "ORIGIN %f %f %f\n", xmin, ymin, zmin);
177 5 : fprintf(fp, "SPACING %f %f %f\n", (xmax - xmin) / (nx - 1), (ymax - ymin) / (ny - 1), dz);
178 :
179 : // Write vector field data
180 5 : fprintf(fp, "\nPOINT_DATA %zu\n", nx * ny * nz);
181 5 : fprintf(fp, "VECTORS %s float\n", field_name);
182 :
183 11 : for (size_t k = 0; k < nz; k++) {
184 57 : for (size_t j = 0; j < ny; j++) {
185 694 : for (size_t i = 0; i < nx; i++) {
186 643 : size_t idx = (k * nx * ny) + IDX_2D(i, j, nx);
187 643 : double w_val = w_data ? w_data[idx] : 0.0;
188 643 : fprintf(fp, "%f %f %f\n", u_data[idx], v_data[idx], w_val);
189 : }
190 : }
191 : }
192 :
193 5 : fclose(fp);
194 : }
195 :
196 4 : void write_vtk_flow_field(const char* filename, const flow_field* field, size_t nx, size_t ny,
197 : size_t nz, double xmin, double xmax, double ymin, double ymax,
198 : double zmin, double zmax) {
199 4 : if (!filename || !field || nx < 2 || ny < 2 || nz < 1 ||
200 2 : xmax <= xmin || ymax <= ymin || (nz > 1 && zmax <= zmin)) {
201 : return;
202 : }
203 :
204 2 : FILE* fp = fopen(filename, "w");
205 2 : if (fp == NULL) {
206 0 : cfd_error("Failed to open VTK flow field output file");
207 0 : return;
208 : }
209 :
210 2 : double dz = (nz > 1) ? (zmax - zmin) / (double)(nz - 1) : 1.0;
211 2 : size_t total = nx * ny * nz;
212 :
213 : // Write VTK header
214 2 : fprintf(fp, "# vtk DataFile Version 3.0\n");
215 2 : fprintf(fp, "CFD Framework Flow Field Output\n");
216 2 : fprintf(fp, "ASCII\n");
217 2 : fprintf(fp, "DATASET STRUCTURED_POINTS\n");
218 2 : fprintf(fp, "DIMENSIONS %zu %zu %zu\n", nx, ny, nz);
219 2 : fprintf(fp, "ORIGIN %f %f %f\n", xmin, ymin, zmin);
220 2 : fprintf(fp, "SPACING %f %f %f\n", (xmax - xmin) / (nx - 1), (ymax - ymin) / (ny - 1), dz);
221 :
222 : // Write multiple field data
223 2 : fprintf(fp, "\nPOINT_DATA %zu\n", total);
224 :
225 : // Velocity vector field
226 2 : fprintf(fp, "VECTORS velocity float\n");
227 4 : for (size_t k = 0; k < nz; k++) {
228 22 : for (size_t j = 0; j < ny; j++) {
229 220 : for (size_t i = 0; i < nx; i++) {
230 200 : size_t idx = (k * nx * ny) + IDX_2D(i, j, nx);
231 200 : double w_val = field->w ? field->w[idx] : 0.0;
232 200 : fprintf(fp, "%f %f %f\n", field->u[idx], field->v[idx], w_val);
233 : }
234 : }
235 : }
236 :
237 : // Pressure scalar field
238 2 : fprintf(fp, "\nSCALARS pressure float 1\n");
239 2 : fprintf(fp, "LOOKUP_TABLE default\n");
240 4 : for (size_t k = 0; k < nz; k++) {
241 22 : for (size_t j = 0; j < ny; j++) {
242 220 : for (size_t i = 0; i < nx; i++) {
243 200 : size_t idx = (k * nx * ny) + IDX_2D(i, j, nx);
244 200 : fprintf(fp, "%f\n", field->p[idx]);
245 : }
246 : }
247 : }
248 :
249 : // Density scalar field
250 2 : fprintf(fp, "\nSCALARS density float 1\n");
251 2 : fprintf(fp, "LOOKUP_TABLE default\n");
252 4 : for (size_t k = 0; k < nz; k++) {
253 22 : for (size_t j = 0; j < ny; j++) {
254 220 : for (size_t i = 0; i < nx; i++) {
255 200 : size_t idx = (k * nx * ny) + IDX_2D(i, j, nx);
256 200 : fprintf(fp, "%f\n", field->rho[idx]);
257 : }
258 : }
259 : }
260 :
261 : // Temperature scalar field
262 2 : fprintf(fp, "\nSCALARS temperature float 1\n");
263 2 : fprintf(fp, "LOOKUP_TABLE default\n");
264 4 : for (size_t k = 0; k < nz; k++) {
265 22 : for (size_t j = 0; j < ny; j++) {
266 220 : for (size_t i = 0; i < nx; i++) {
267 200 : size_t idx = (k * nx * ny) + IDX_2D(i, j, nx);
268 200 : fprintf(fp, "%f\n", field->T[idx]);
269 : }
270 : }
271 : }
272 :
273 2 : fclose(fp);
274 : }
275 :
276 : // New run-based output functions
277 0 : static void get_run_filepath(char* buffer, size_t buffer_size, const char* filename) {
278 0 : char run_dir[512];
279 0 : cfd_get_run_directory(run_dir, sizeof(run_dir));
280 :
281 : // If no run directory exists yet, create one
282 0 : if (strlen(run_dir) == 0) {
283 0 : cfd_create_run_directory(run_dir, sizeof(run_dir));
284 : }
285 :
286 : // Build full path
287 : #ifdef _WIN32
288 : snprintf(buffer, buffer_size, "%s\\%s", run_dir, filename);
289 : #else
290 0 : snprintf(buffer, buffer_size, "%s/%s", run_dir, filename);
291 : #endif
292 0 : }
293 :
294 0 : void write_vtk_output_run(const char* filename, const char* field_name, const double* data,
295 : size_t nx, size_t ny, size_t nz, double xmin, double xmax, double ymin,
296 : double ymax, double zmin, double zmax) {
297 0 : char filepath[1024];
298 0 : get_run_filepath(filepath, sizeof(filepath), filename);
299 0 : write_vtk_output(filepath, field_name, data, nx, ny, nz, xmin, xmax, ymin, ymax, zmin, zmax);
300 0 : }
301 :
302 0 : void write_vtk_vector_output_run(const char* filename, const char* field_name, const double* u_data,
303 : const double* v_data, const double* w_data, size_t nx, size_t ny,
304 : size_t nz, double xmin, double xmax, double ymin, double ymax,
305 : double zmin, double zmax) {
306 0 : char filepath[1024];
307 0 : get_run_filepath(filepath, sizeof(filepath), filename);
308 0 : write_vtk_vector_output(filepath, field_name, u_data, v_data, w_data, nx, ny, nz, xmin, xmax,
309 : ymin, ymax, zmin, zmax);
310 0 : }
311 :
312 0 : void write_vtk_flow_field_run(const char* filename, const flow_field* field, size_t nx, size_t ny,
313 : size_t nz, double xmin, double xmax, double ymin, double ymax,
314 : double zmin, double zmax) {
315 0 : char filepath[1024];
316 0 : get_run_filepath(filepath, sizeof(filepath), filename);
317 0 : write_vtk_flow_field(filepath, field, nx, ny, nz, xmin, xmax, ymin, ymax, zmin, zmax);
318 0 : }
|