FrontISTR  5.9.0
Large-scale structural analysis program with finit element method
hecmw_fstr_output_exo.c
Go to the documentation of this file.
1 /*****************************************************************************
2  * Copyright (c) 2019 FrontISTR Commons
3  * This software is released under the MIT License, see LICENSE.txt
4  *
5  * Exodus II output via NetCDF API (without libexodus dependency)
6  *****************************************************************************/
8 
9 #ifdef WITH_NETCDF
10 
11 #include <stdio.h>
12 #include <stdlib.h>
13 #include <string.h>
14 #include <time.h>
15 #include <netcdf.h>
16 #include "hecmw_malloc.h"
17 #include "hecmw_etype.h"
18 #include "hecmw_vis_mem_util.h"
19 #include "hecmw_vis_comm_util.h"
20 
21 /*============================================================================
22  * Constants
23  *============================================================================*/
24 #define EXO_MAX_STR_LENGTH 33
25 #define EXO_MAX_LINE_LENGTH 81
26 #define EXO_API_VERSION 8.03f
27 #define EXO_DB_VERSION 8.03f
28 #define EXO_FILE_SIZE 1
29 
30 /* Maximum number of distinct element types we handle */
31 #define MAX_ELEM_BLOCKS 64
32 
33 /*============================================================================
34  * FrontISTR element type -> Exodus element type string mapping
35  *
36  * FrontISTR uses numeric element type codes defined in hecmw_common_define.h.
37  * Exodus requires a string like "TETRA4".
38  *
39  * The actual number of nodes written to connectivity is determined at runtime
40  * from elem_node_index (with shift applied), matching the VTK output behavior
41  * exactly. The exo_name here tells ParaView how to interpret the topology.
42  *
43  * Element categories:
44  * 1xx: Rod/Bar 2xx: 2D surface 3xx: 3D solid
45  * 4xx: Master/slave 5xx: Joint 6xx: Beam
46  * 7xx: Shell 9xx: Interface(LN) 10xx: Patch(PT)
47  *============================================================================*/
48 typedef struct {
49  int hecmw_type; /* FrontISTR element type code */
50  const char *exo_name; /* Exodus element type string */
51 } ElemTypeMap;
52 
53 static const ElemTypeMap elem_type_table[] = {
54  /* --- 1D Rod/Bar elements --- */
55  {111, "BAR2"}, /* ROD1: 2 nodes */
56  {112, "BAR3"}, /* ROD2: 3 nodes */
57  {301, "BAR2"}, /* ROD31: 2 nodes */
58  /* --- 2D Surface elements --- */
59  {231, "TRI3"}, /* TRI1: 3 nodes */
60  {232, "TRI6"}, /* TRI2: 6 nodes */
61  {2322, "TRI3"}, /* TRI22: variant, VTK treats as TRI */
62  {241, "QUAD4"}, /* QUA1: 4 nodes */
63  {242, "QUAD8"}, /* QUA2: 8 nodes */
64  /* --- 3D Solid elements --- */
65  {341, "TETRA4"}, /* TET1: 4 nodes */
66  {3414, "TETRA4"}, /* TET1_4: 4 nodes (variant) */
67  {342, "TETRA10"}, /* TET2: 10 nodes (node reorder needed) */
68  {3422, "TETRA4"}, /* TET22: VTK treats as TETRA (4 nodes) */
69  {351, "WEDGE6"}, /* PRI1: 6 nodes */
70  {352, "WEDGE15"}, /* PRI2: 15 nodes */
71  {361, "HEX8"}, /* HEX1: 8 nodes */
72  {3614, "HEX8"}, /* HEX1_4: 8 nodes (variant) */
73  {362, "HEX20"}, /* HEX2: 20 nodes */
74  {371, "PYRAMID5"}, /* PYR1: 5 nodes */
75  {372, "PYRAMID13"}, /* PYR2: 13 nodes */
76  /* --- Master/Slave surface elements --- */
77  {431, "TETRA4"}, /* MST1: 4 nodes, VTK=TRI but stored as tet */
78  {432, "TETRA4"}, /* MST2: 7 nodes, VTK=TRI */
79  {441, "PYRAMID5"}, /* MSQ1: 5 nodes, VTK=QUAD but stored as pyr */
80  {442, "PYRAMID5"}, /* MSQ2: 9 nodes, VTK=QUAD */
81  /* --- Joint elements --- */
82  {501, "BAR2"}, /* JTB1: 2 nodes */
83  {511, "BAR2"}, /* SPGDPT1: 2 nodes */
84  {531, "TRI3"}, /* JTT1: 6 nodes total, VTK=TRI (uses first 3) */
85  {532, "TRI3"}, /* JTT2: 12 nodes total, VTK=TRI */
86  {541, "QUAD4"}, /* JTQ1: 8 nodes total, VTK=QUAD (uses first 4) */
87  {542, "QUAD4"}, /* JTQ2: 16 nodes total, VTK=QUAD */
88  /* --- Beam elements --- */
89  {611, "BAR2"}, /* BEM1: 2 nodes */
90  {612, "BAR3"}, /* BEM2: 3 nodes */
91  {641, "BAR2"}, /* BEM3: 4 nodes, shift=2 -> 2 nodes (mixed beam-341) */
92  /* --- Shell elements --- */
93  {731, "TRI3"}, /* SHT1: 3 nodes */
94  {732, "TRI6"}, /* SHT2: 6 nodes */
95  {741, "QUAD4"}, /* SHQ1: 4 nodes */
96  {742, "QUAD8"}, /* SHQ2: 8 nodes */
97  {743, "QUAD9"}, /* SHQ3: 9 nodes */
98  {761, "TRI3"}, /* SHT6: 6 nodes, shift=3 -> 3 nodes (mixed shell-solid) */
99  {781, "QUAD4"}, /* SHQ8: 8 nodes, shift=4 -> 4 nodes (mixed shell-solid) */
100  /* --- Interface/Link elements (LN) --- all 2-node lines */
101  {911, "BAR2"}, {912, "BAR2"}, {913, "BAR2"}, {914, "BAR2"}, {915, "BAR2"}, {916, "BAR2"},
102  {921, "BAR2"}, {922, "BAR2"}, {923, "BAR2"}, {924, "BAR2"}, {925, "BAR2"}, {926, "BAR2"},
103  {931, "BAR2"}, {932, "BAR2"}, {933, "BAR2"}, {934, "BAR2"}, {935, "BAR2"}, {936, "BAR2"},
104  {941, "BAR2"}, {942, "BAR2"}, {943, "BAR2"}, {944, "BAR2"}, {945, "BAR2"}, {946, "BAR2"},
105  {951, "BAR2"}, {952, "BAR2"}, {953, "BAR2"}, {954, "BAR2"}, {955, "BAR2"}, {956, "BAR2"},
106  {961, "BAR2"}, {962, "BAR2"}, {963, "BAR2"}, {964, "BAR2"}, {965, "BAR2"}, {966, "BAR2"},
107  /* --- Patch elements (PT) --- */
108  {1031, "TRI3"}, /* PTT1: 3 nodes */
109  {1032, "TRI6"}, /* PTT2: 6 nodes */
110  {1041, "QUAD4"}, /* PTQ1: 4 nodes */
111  {1042, "QUAD8"}, /* PTQ2: 8 nodes */
112  /* --- Sentinel --- */
113  {0, NULL}
114 };
115 
116 /*============================================================================
117  * Helper: NetCDF error check macro
118  *============================================================================*/
119 #define NCERR(e) do { \
120  if ((e) != NC_NOERR) { \
121  fprintf(stderr, "NetCDF error at %s:%d: %s\n", \
122  __FILE__, __LINE__, nc_strerror(e)); \
123  return; \
124  } \
125 } while(0)
126 
127 /*============================================================================
128  * Helper: look up Exodus element type string for a FrontISTR type code
129  *============================================================================*/
130 static const ElemTypeMap* get_elem_type_map(int hecmw_type)
131 {
132  int i;
133  for (i = 0; elem_type_table[i].exo_name != NULL; i++) {
134  if (elem_type_table[i].hecmw_type == hecmw_type)
135  return &elem_type_table[i];
136  }
137  return NULL;
138 }
139 
140 /*============================================================================
141  * Helper: get node shift for special element types (same as VTK output)
142  *============================================================================*/
143 static long long get_node_shift(int hecmw_type)
144 {
145  switch (hecmw_type) {
146  case 641: return 2;
147  case 761: return 3;
148  case 781: return 4;
149  default: return 0;
150  }
151 }
152 
153 /*============================================================================
154  * Helper: get the effective number of nodes for an element
155  * (after applying shift, matching VTK output behavior)
156  *============================================================================*/
157 static int get_effective_n_nodes(struct hecmwST_local_mesh *mesh, int elem_idx)
158 {
159  long long jS = mesh->elem_node_index[elem_idx];
160  long long jE = mesh->elem_node_index[elem_idx + 1];
161  long long shift = get_node_shift(mesh->elem_type[elem_idx]);
162  return jE - jS - shift;
163 }
164 
165 /*============================================================================
166  * Element block info (gathered by scanning mesh element types)
167  *============================================================================*/
168 typedef struct {
169  int hecmw_type; /* FrontISTR element type code */
170  const char *exo_name; /* Exodus element type string */
171  int nod_per_elem; /* nodes per element */
172  int num_elem; /* number of elements in this block */
173  int *elem_indices; /* indices into original element array */
174 } ElemBlockInfo;
175 
176 /*============================================================================
177  * Helper: write a padded string into a char array for NetCDF
178  *============================================================================*/
179 static void pad_string(char *buf, int len, const char *str)
180 {
181  int slen = (int)strlen(str);
182  int i;
183  if (slen > len) slen = len;
184  for (i = 0; i < slen; i++) buf[i] = str[i];
185  for (i = slen; i < len; i++) buf[i] = '\0';
186 }
187 
188 /*============================================================================
189  * Helper: expand multi-component FrontISTR variable labels into per-component
190  * Exodus variable names.
191  *
192  * FrontISTR stores a single label like "DISPLACEMENT" with nn_dof[i]=3,
193  * meaning 3 components (X,Y,Z). Exodus needs individual names per component:
194  * "DISPLACEMENT_X", "DISPLACEMENT_Y", "DISPLACEMENT_Z"
195  *
196  * For 6-component tensor data (e.g., stress):
197  * "STRESS_XX", "STRESS_YY", "STRESS_ZZ", "STRESS_XY", "STRESS_YZ", "STRESS_ZX"
198  *============================================================================*/
199 static const char *comp_suffix_1[] = {""};
200 static const char *comp_suffix_2[] = {"_X", "_Y"};
201 static const char *comp_suffix_3[] = {"_X", "_Y", "_Z"};
202 static const char *comp_suffix_3_principal[] = {"_1st", "_2nd", "_3rd"};
203 static const char *comp_suffix_6[] = {"_XX", "_YY", "_ZZ", "_XY", "_YZ", "_ZX"};
204 static const char *comp_suffix_7[] = {"_XX", "_YY", "_ZZ", "_XY", "_YZ", "_ZX", "_MISES"};
205 
206 static int label_is_principal(const char *label)
207 {
208  /* Match FrontISTR labels:
209  * NodalPrincipalSTRESS, NodalPrincipalSTRAIN,
210  * ElementalPrincipalSTRESS, ElementalPrincipalSTRAIN, etc. */
211  return (strstr(label, "Principal") != NULL);
212 }
213 
214 static const char** get_comp_suffixes(int ndof, const char *label)
215 {
216  switch (ndof) {
217  case 1: return comp_suffix_1;
218  case 2: return comp_suffix_2;
219  case 3: return label_is_principal(label) ? comp_suffix_3_principal
220  : comp_suffix_3;
221  case 6: return comp_suffix_6;
222  case 7: return comp_suffix_7;
223  default: return NULL;
224  }
225 }
226 
227 /*============================================================================
228  * Main Exodus output function
229  *
230  * This function writes an Exodus II file using NetCDF API calls directly.
231  * It follows the same calling convention as vtk_output / bin_vtk_output.
232  *
233  * On the first call (time_step_count == 0), it creates the file and defines
234  * the full schema. On subsequent calls, it opens the file and appends a
235  * new time step.
236  *
237  * File naming:
238  * - Each rank writes: {outfile1}/{outfile}.{myrank}.exo
239  * (No pvtu-equivalent needed: ParaView can load .exo.N.M Nemesis files,
240  * or individual per-rank .exo files can be loaded as a group.)
241  *============================================================================*/
242 void exodus_output(struct hecmwST_local_mesh *mesh,
243  struct hecmwST_result_data *data,
244  char *outfile, char *outfile1,
245  HECMW_Comm VIS_COMM,
246  int per_step)
247 {
248  int i, j, k;
249  long long jS, jE;
250  int myrank, petot;
251  int n_node, n_elem, shift;
252  int data_tot_n, data_tot_e;
253  int table342[10] = {0, 1, 2, 3, 6, 4, 5, 7, 8, 9};
254  char *p;
255 
256  int ncid, stat;
257  int old_fill_mode;
258  static int time_step_count = 0;
259 
260  /*
261  * Two modes:
262  * per_step == 0 (EXODUS): All timesteps in one file (append mode)
263  * per_step == 1 (STEP_EXODUS): One file per timestep (always create new)
264  *
265  * For per_step==0, the file path is fixed across calls (static).
266  * For per_step==1, a new file is created each call with the step number.
267  */
268  static char file_exo[HECMW_FILENAME_LEN];
269  static char dir_exo[HECMW_FILENAME_LEN];
270  static char basename_exo[HECMW_FILENAME_LEN];
271 
272  /* --- Element block working data --- */
273  int num_blocks = 0;
274  ElemBlockInfo blocks[MAX_ELEM_BLOCKS];
275 
276  /* --- Count total DOFs for nodal/element data --- */
277  int total_nod_var = 0; /* total number of scalar Exodus nodal variables */
278  int total_elem_var = 0; /* total number of scalar Exodus element variables */
279 
280  HECMW_Comm_rank(VIS_COMM, &myrank);
281  HECMW_Comm_size(VIS_COMM, &petot);
282  n_node = mesh->n_node;
283  n_elem = mesh->n_elem;
284 
285  data_tot_n = 0;
286  for (i = 0; i < data->nn_component; i++)
287  data_tot_n += data->nn_dof[i];
288 
289  data_tot_e = 0;
290  for (i = 0; i < data->ne_component; i++)
291  data_tot_e += data->ne_dof[i];
292 
293  /* Count total scalar variables for Exodus (expand multi-component) */
294  for (i = 0; i < data->nn_component; i++)
295  total_nod_var += data->nn_dof[i];
296  for (i = 0; i < data->ne_component; i++)
297  total_elem_var += data->ne_dof[i];
298 
299  /* --- Build file path --- */
300  if (time_step_count == 0) {
301  /*
302  * First call: compute dir_exo and basename_exo from caller's args.
303  *
304  * outfile1 is e.g. "vis_out_psf.0001" (with timestep suffix).
305  * We strip the suffix to get a stable directory: "vis_out_psf_exo".
306  */
307 
308  /* Build directory: strip timestep suffix from outfile1 and append _exo */
309  snprintf(dir_exo, sizeof(dir_exo), "%s", outfile1);
310  p = strrchr(dir_exo, '.');
311  if (p != NULL) *p = '\0'; /* "vis_out_psf.0001" -> "vis_out_psf" */
312  strncat(dir_exo, "_exo", sizeof(dir_exo) - strlen(dir_exo) - 1);
313  /* dir_exo = "vis_out_psf_exo" */
314 
315  /* Build basename: strip directory prefix and timestep suffix from outfile */
316  snprintf(basename_exo, sizeof(basename_exo), "%s", outfile);
317  p = strrchr(basename_exo, '/');
318  if (p != NULL) {
319  memmove(basename_exo, p + 1, strlen(p + 1) + 1);
320  }
321  p = strrchr(basename_exo, '.');
322  if (p != NULL) *p = '\0'; /* "vis_out_psf.0001" -> "vis_out_psf" */
323  }
324 
325  if (per_step) {
326  /*
327  * STEP_EXODUS mode: one file per timestep
328  * serial: {basename}.step{NNNN}.exo
329  * parallel: {basename}.step{NNNN}.exo.{nprocs}.{rank} (Nemesis)
330  */
331  if (petot == 1) {
332  snprintf(file_exo, sizeof(file_exo), "%s/%s.step%04d.exo",
333  dir_exo, basename_exo, time_step_count);
334  } else {
335  snprintf(file_exo, sizeof(file_exo), "%s/%s.step%04d.exo.%d.%d",
336  dir_exo, basename_exo, time_step_count, petot, myrank);
337  }
338  } else {
339  /*
340  * EXODUS mode: all timesteps in one file (only set path on first call)
341  * serial: {basename}.exo
342  * parallel: {basename}.exo.{nprocs}.{rank} (Nemesis)
343  */
344  if (time_step_count == 0) {
345  if (petot == 1) {
346  snprintf(file_exo, sizeof(file_exo), "%s/%s.exo",
347  dir_exo, basename_exo);
348  } else {
349  snprintf(file_exo, sizeof(file_exo), "%s/%s.exo.%d.%d",
350  dir_exo, basename_exo, petot, myrank);
351  }
352  }
353  /* For subsequent calls, file_exo is already set (static) */
354  }
355 
356  /* Create the directory (needed on first call, or every call for per_step) */
357  if (time_step_count == 0 || per_step) {
358  if (HECMW_ctrl_make_subdir(file_exo)) {
359  HECMW_vis_print_exit("ERROR: HEC-MW-VIS-E0009: Cannot open output directory for Exodus");
360  }
361  }
362 
363  /*========================================================================
364  * PASS 1: Scan elements and group into blocks by element type
365  *========================================================================*/
366  memset(blocks, 0, sizeof(blocks));
367  for (i = 0; i < n_elem; i++) {
368  int etype = mesh->elem_type[i];
369  int found = 0;
370  for (j = 0; j < num_blocks; j++) {
371  if (blocks[j].hecmw_type == etype) {
372  blocks[j].num_elem++;
373  found = 1;
374  break;
375  }
376  }
377  if (!found) {
378  if (num_blocks >= MAX_ELEM_BLOCKS) {
379  HECMW_vis_print_exit("ERROR: Too many element types for Exodus output");
380  }
381  const ElemTypeMap *map = get_elem_type_map(etype);
382  if (map == NULL) {
383  fprintf(stderr, "WARNING: Unknown element type %d, skipping\n", etype);
384  continue;
385  }
386  blocks[num_blocks].hecmw_type = etype;
387  blocks[num_blocks].exo_name = map->exo_name;
388  /* Determine nodes per element from first element of this type */
389  blocks[num_blocks].nod_per_elem = get_effective_n_nodes(mesh, i);
390  blocks[num_blocks].num_elem = 1;
391  num_blocks++;
392  }
393  }
394 
395  /* Allocate index arrays and fill element indices for each block */
396  for (j = 0; j < num_blocks; j++) {
397  blocks[j].elem_indices = (int *)HECMW_malloc(sizeof(int) * blocks[j].num_elem);
398  if (blocks[j].elem_indices == NULL)
399  HECMW_vis_print_exit("HECMW_malloc failed for block elem_indices");
400  blocks[j].num_elem = 0; /* reset count, will re-fill */
401  }
402  for (i = 0; i < n_elem; i++) {
403  int etype = mesh->elem_type[i];
404  for (j = 0; j < num_blocks; j++) {
405  if (blocks[j].hecmw_type == etype) {
406  blocks[j].elem_indices[blocks[j].num_elem] = i;
407  blocks[j].num_elem++;
408  break;
409  }
410  }
411  }
412 
413  /*========================================================================
414  * CREATE or OPEN file
415  *
416  * per_step==1: Always create a new file (one file per timestep)
417  * per_step==0: Create on first call, open for append on subsequent calls
418  *========================================================================*/
419  if (per_step || time_step_count == 0) {
420  /*====================================================================
421  * FIRST CALL: Create file and define full schema
422  *====================================================================*/
423  stat = nc_create(file_exo, NC_CLOBBER | NC_64BIT_OFFSET, &ncid);
424  NCERR(stat);
425  stat = nc_set_fill(ncid, NC_NOFILL, &old_fill_mode);
426  NCERR(stat);
427 
428  /*----------------------------------------------------------------
429  * Global Attributes [REQUIRED]
430  *----------------------------------------------------------------*/
431  {
432  float api_ver = EXO_API_VERSION;
433  float db_ver = EXO_DB_VERSION;
434  int word_sz = (int)sizeof(double);
435  int file_sz = EXO_FILE_SIZE;
436  char title[] = "FrontISTR Exodus output";
437 
438  stat = nc_put_att_float(ncid, NC_GLOBAL, "api_version", NC_FLOAT, 1, &api_ver);
439  NCERR(stat);
440  stat = nc_put_att_float(ncid, NC_GLOBAL, "version", NC_FLOAT, 1, &db_ver);
441  NCERR(stat);
442  stat = nc_put_att_int(ncid, NC_GLOBAL, "floating_point_word_size", NC_INT, 1, &word_sz);
443  NCERR(stat);
444  stat = nc_put_att_int(ncid, NC_GLOBAL, "file_size", NC_INT, 1, &file_sz);
445  NCERR(stat);
446  stat = nc_put_att_text(ncid, NC_GLOBAL, "title", strlen(title), title);
447  NCERR(stat);
448  }
449 
450  /*----------------------------------------------------------------
451  * Dimensions [REQUIRED]
452  *----------------------------------------------------------------*/
453  int dim_len_string, dim_len_name, dim_len_line, dim_four;
454  int dim_time_step, dim_num_dim, dim_num_nodes, dim_num_elem;
455  int dim_num_el_blk;
456  int dim_num_qa_rec;
457 
458  stat = nc_def_dim(ncid, "len_string", EXO_MAX_STR_LENGTH, &dim_len_string);
459  NCERR(stat);
460  stat = nc_def_dim(ncid, "len_name", EXO_MAX_STR_LENGTH, &dim_len_name);
461  NCERR(stat);
462  stat = nc_def_dim(ncid, "len_line", EXO_MAX_LINE_LENGTH, &dim_len_line);
463  NCERR(stat);
464  stat = nc_def_dim(ncid, "four", 4, &dim_four);
465  NCERR(stat);
466  stat = nc_def_dim(ncid, "time_step", NC_UNLIMITED, &dim_time_step);
467  NCERR(stat);
468  stat = nc_def_dim(ncid, "num_dim", 3, &dim_num_dim);
469  NCERR(stat);
470  stat = nc_def_dim(ncid, "num_nodes", (size_t)n_node, &dim_num_nodes);
471  NCERR(stat);
472  stat = nc_def_dim(ncid, "num_elem", (size_t)n_elem, &dim_num_elem);
473  NCERR(stat);
474  stat = nc_def_dim(ncid, "num_el_blk", (size_t)num_blocks, &dim_num_el_blk);
475  NCERR(stat);
476 
477  /* [OPTIONAL] QA record dimension */
478  stat = nc_def_dim(ncid, "num_qa_rec", 1, &dim_num_qa_rec);
479  NCERR(stat);
480 
481  /* Per-block dimensions */
482  int dim_el_in_blk[MAX_ELEM_BLOCKS];
483  int dim_nod_per_el[MAX_ELEM_BLOCKS];
484  for (j = 0; j < num_blocks; j++) {
485  char dname[EXO_MAX_STR_LENGTH];
486  snprintf(dname, sizeof(dname), "num_el_in_blk%d", j + 1);
487  stat = nc_def_dim(ncid, dname, (size_t)blocks[j].num_elem, &dim_el_in_blk[j]);
488  NCERR(stat);
489  snprintf(dname, sizeof(dname), "num_nod_per_el%d", j + 1);
490  stat = nc_def_dim(ncid, dname, (size_t)blocks[j].nod_per_elem, &dim_nod_per_el[j]);
491  NCERR(stat);
492  }
493 
494  /* Variable count dimensions */
495  int dim_num_nod_var = -1, dim_num_elem_var = -1;
496  if (total_nod_var > 0) {
497  stat = nc_def_dim(ncid, "num_nod_var", (size_t)total_nod_var, &dim_num_nod_var);
498  NCERR(stat);
499  }
500  if (total_elem_var > 0) {
501  stat = nc_def_dim(ncid, "num_elem_var", (size_t)total_elem_var, &dim_num_elem_var);
502  NCERR(stat);
503  }
504 
505  /*----------------------------------------------------------------
506  * Variables: Time [REQUIRED]
507  *----------------------------------------------------------------*/
508  int var_time_whole;
509  stat = nc_def_var(ncid, "time_whole", NC_DOUBLE, 1, &dim_time_step, &var_time_whole);
510  NCERR(stat);
511 
512  /*----------------------------------------------------------------
513  * Variables: Coordinates [REQUIRED]
514  *----------------------------------------------------------------*/
515  int var_coordx, var_coordy, var_coordz;
516  stat = nc_def_var(ncid, "coordx", NC_DOUBLE, 1, &dim_num_nodes, &var_coordx);
517  NCERR(stat);
518  stat = nc_def_var(ncid, "coordy", NC_DOUBLE, 1, &dim_num_nodes, &var_coordy);
519  NCERR(stat);
520  stat = nc_def_var(ncid, "coordz", NC_DOUBLE, 1, &dim_num_nodes, &var_coordz);
521  NCERR(stat);
522 
523  /* coor_names [REQUIRED] */
524  int var_coor_names;
525  {
526  int dims2[2] = {dim_num_dim, dim_len_name};
527  stat = nc_def_var(ncid, "coor_names", NC_CHAR, 2, dims2, &var_coor_names);
528  NCERR(stat);
529  }
530 
531  /*----------------------------------------------------------------
532  * Variables: Element Block info [REQUIRED]
533  *----------------------------------------------------------------*/
534  int var_eb_status, var_eb_prop1;
535  stat = nc_def_var(ncid, "eb_status", NC_INT, 1, &dim_num_el_blk, &var_eb_status);
536  NCERR(stat);
537  stat = nc_def_var(ncid, "eb_prop1", NC_INT, 1, &dim_num_el_blk, &var_eb_prop1);
538  NCERR(stat);
539  stat = nc_put_att_text(ncid, var_eb_prop1, "name", 2, "ID");
540  NCERR(stat);
541 
542  /* [OPTIONAL] eb_names */
543  int var_eb_names;
544  {
545  int dims2[2] = {dim_num_el_blk, dim_len_name};
546  stat = nc_def_var(ncid, "eb_names", NC_CHAR, 2, dims2, &var_eb_names);
547  NCERR(stat);
548  }
549 
550  /* Connectivity per block [REQUIRED] */
551  int var_connect[MAX_ELEM_BLOCKS];
552  for (j = 0; j < num_blocks; j++) {
553  char vname[EXO_MAX_STR_LENGTH];
554  snprintf(vname, sizeof(vname), "connect%d", j + 1);
555  int dims2[2] = {dim_el_in_blk[j], dim_nod_per_el[j]};
556  stat = nc_def_var(ncid, vname, NC_INT, 2, dims2, &var_connect[j]);
557  NCERR(stat);
558  stat = nc_put_att_text(ncid, var_connect[j], "elem_type",
559  strlen(blocks[j].exo_name), blocks[j].exo_name);
560  NCERR(stat);
561  }
562 
563  /*----------------------------------------------------------------
564  * Variables: QA records [OPTIONAL - nice to have]
565  *----------------------------------------------------------------*/
566  int var_qa_records;
567  {
568  int dims3[3] = {dim_num_qa_rec, dim_four, dim_len_string};
569  stat = nc_def_var(ncid, "qa_records", NC_CHAR, 3, dims3, &var_qa_records);
570  NCERR(stat);
571  }
572 
573  /*----------------------------------------------------------------
574  * Variables: ID maps [OPTIONAL - nice to have for parallel]
575  *----------------------------------------------------------------*/
576  int var_node_num_map, var_elem_num_map;
577  stat = nc_def_var(ncid, "node_num_map", NC_INT, 1, &dim_num_nodes, &var_node_num_map);
578  NCERR(stat);
579  stat = nc_def_var(ncid, "elem_num_map", NC_INT, 1, &dim_num_elem, &var_elem_num_map);
580  NCERR(stat);
581 
582  /*----------------------------------------------------------------
583  * Variables: Nodal result variable names and data [REQUIRED]
584  *----------------------------------------------------------------*/
585  int var_name_nod_var = -1;
586  int var_vals_nod[512]; /* enough for many variables */
587  if (total_nod_var > 0) {
588  int dims2[2] = {dim_num_nod_var, dim_len_name};
589  stat = nc_def_var(ncid, "name_nod_var", NC_CHAR, 2, dims2, &var_name_nod_var);
590  NCERR(stat);
591 
592  for (i = 0; i < total_nod_var; i++) {
593  char vname[EXO_MAX_STR_LENGTH];
594  snprintf(vname, sizeof(vname), "vals_nod_var%d", i + 1);
595  int dims_ts_nn[2] = {dim_time_step, dim_num_nodes};
596  stat = nc_def_var(ncid, vname, NC_DOUBLE, 2, dims_ts_nn, &var_vals_nod[i]);
597  NCERR(stat);
598  }
599  }
600 
601  /*----------------------------------------------------------------
602  * Variables: Element result variable names and data [REQUIRED]
603  *----------------------------------------------------------------*/
604  int var_name_elem_var = -1;
605  int var_vals_elem[512]; /* [var_idx * num_blocks + blk_idx] */
606  int var_elem_var_tab = -1;
607  if (total_elem_var > 0) {
608  int dims2[2] = {dim_num_elem_var, dim_len_name};
609  stat = nc_def_var(ncid, "name_elem_var", NC_CHAR, 2, dims2, &var_name_elem_var);
610  NCERR(stat);
611 
612  /* Truth table */
613  {
614  int dims_tt[2] = {dim_num_el_blk, dim_num_elem_var};
615  stat = nc_def_var(ncid, "elem_var_tab", NC_INT, 2, dims_tt, &var_elem_var_tab);
616  NCERR(stat);
617  }
618 
619  /* vals_elem_var{v}eb{b} for each variable and block */
620  for (i = 0; i < total_elem_var; i++) {
621  for (j = 0; j < num_blocks; j++) {
622  char vname[EXO_MAX_STR_LENGTH];
623  snprintf(vname, sizeof(vname), "vals_elem_var%deb%d", i + 1, j + 1);
624  int dims_ts_eb[2] = {dim_time_step, dim_el_in_blk[j]};
625  stat = nc_def_var(ncid, vname, NC_DOUBLE, 2, dims_ts_eb,
626  &var_vals_elem[i * num_blocks + j]);
627  NCERR(stat);
628  }
629  }
630  }
631 
632  /*================================================================
633  * END DEFINE MODE
634  *================================================================*/
635  stat = nc_enddef(ncid);
636  NCERR(stat);
637 
638  /*================================================================
639  * WRITE STATIC DATA (coordinates, connectivity, maps, etc.)
640  * These are written only once.
641  *================================================================*/
642 
643  /* --- Coordinates (XYZ separated) --- */
644  {
645  double *cx = (double *)HECMW_malloc(sizeof(double) * n_node);
646  double *cy = (double *)HECMW_malloc(sizeof(double) * n_node);
647  double *cz = (double *)HECMW_malloc(sizeof(double) * n_node);
648  if (cx == NULL || cy == NULL || cz == NULL) {
649  HECMW_free(cx); HECMW_free(cy); HECMW_free(cz);
650  HECMW_vis_print_exit("HECMW_malloc failed for coordinate arrays");
651  }
652  for (i = 0; i < n_node; i++) {
653  cx[i] = mesh->node[3 * i];
654  cy[i] = mesh->node[3 * i + 1];
655  cz[i] = mesh->node[3 * i + 2];
656  }
657  stat = nc_put_var_double(ncid, var_coordx, cx); NCERR(stat);
658  stat = nc_put_var_double(ncid, var_coordy, cy); NCERR(stat);
659  stat = nc_put_var_double(ncid, var_coordz, cz); NCERR(stat);
660  HECMW_free(cx); HECMW_free(cy); HECMW_free(cz);
661  }
662 
663  /* --- Coordinate names --- */
664  {
665  char names[3][EXO_MAX_STR_LENGTH];
666  pad_string(names[0], EXO_MAX_STR_LENGTH, "x");
667  pad_string(names[1], EXO_MAX_STR_LENGTH, "y");
668  pad_string(names[2], EXO_MAX_STR_LENGTH, "z");
669  stat = nc_put_var_text(ncid, var_coor_names, &names[0][0]);
670  NCERR(stat);
671  }
672 
673  /* --- Element block status and IDs --- */
674  {
675  int *eb_stat = (int *)HECMW_malloc(sizeof(int) * num_blocks);
676  int *eb_ids = (int *)HECMW_malloc(sizeof(int) * num_blocks);
677  if (eb_stat == NULL || eb_ids == NULL) {
678  HECMW_free(eb_stat); HECMW_free(eb_ids);
679  HECMW_vis_print_exit("HECMW_malloc failed for element block arrays");
680  }
681  for (j = 0; j < num_blocks; j++) {
682  eb_stat[j] = 1; /* active */
683  eb_ids[j] = j + 1; /* block ID = 1,2,3,... */
684  }
685  stat = nc_put_var_int(ncid, var_eb_status, eb_stat); NCERR(stat);
686  stat = nc_put_var_int(ncid, var_eb_prop1, eb_ids); NCERR(stat);
687  HECMW_free(eb_stat); HECMW_free(eb_ids);
688  }
689 
690  /* --- [OPTIONAL] Element block names --- */
691  {
692  char *names_buf = (char *)HECMW_malloc(num_blocks * EXO_MAX_STR_LENGTH);
693  if (names_buf == NULL)
694  HECMW_vis_print_exit("HECMW_malloc failed for block names");
695  for (j = 0; j < num_blocks; j++) {
696  char label[EXO_MAX_STR_LENGTH];
697  snprintf(label, sizeof(label), "block_%d_%s", j + 1, blocks[j].exo_name);
698  pad_string(&names_buf[j * EXO_MAX_STR_LENGTH], EXO_MAX_STR_LENGTH, label);
699  }
700  stat = nc_put_var_text(ncid, var_eb_names, names_buf);
701  NCERR(stat);
702  HECMW_free(names_buf);
703  }
704 
705  /* --- Connectivity per block --- */
706  for (j = 0; j < num_blocks; j++) {
707  int ne = blocks[j].num_elem;
708  int nn = blocks[j].nod_per_elem;
709  int etype = blocks[j].hecmw_type;
710  int node_shift = get_node_shift(etype);
711  int *conn = (int *)HECMW_malloc(sizeof(int) * ne * nn);
712  if (conn == NULL)
713  HECMW_vis_print_exit("HECMW_malloc failed for connectivity");
714 
715  for (i = 0; i < ne; i++) {
716  int ei = blocks[j].elem_indices[i]; /* original element index */
717  jS = mesh->elem_node_index[ei];
718  jE = mesh->elem_node_index[ei + 1];
719 
720  if (etype == 342) {
721  /* TET10: reorder nodes (same table as VTK output) */
722  for (k = 0; k < nn; k++) {
723  conn[i * nn + k] = mesh->elem_node_item[jS + table342[k]];
724  /* NOTE: Exodus uses 1-based indexing (no -1) */
725  }
726  } else {
727  for (k = 0; k < nn; k++) {
728  conn[i * nn + k] = mesh->elem_node_item[jS + k];
729  }
730  }
731  }
732  stat = nc_put_var_int(ncid, var_connect[j], conn);
733  NCERR(stat);
734  HECMW_free(conn);
735  }
736 
737  /* --- [OPTIONAL] QA records --- */
738  {
739  char qa_data[1][4][EXO_MAX_STR_LENGTH];
740  time_t now = time(NULL);
741  struct tm *t = localtime(&now);
742  char datestr[32], timestr[32];
743  strftime(datestr, sizeof(datestr), "%Y%m%d", t);
744  strftime(timestr, sizeof(timestr), "%H:%M:%S", t);
745 
746  pad_string(qa_data[0][0], EXO_MAX_STR_LENGTH, "FrontISTR");
747  pad_string(qa_data[0][1], EXO_MAX_STR_LENGTH, "1.0");
748  pad_string(qa_data[0][2], EXO_MAX_STR_LENGTH, datestr);
749  pad_string(qa_data[0][3], EXO_MAX_STR_LENGTH, timestr);
750  stat = nc_put_var_text(ncid, var_qa_records, &qa_data[0][0][0]);
751  NCERR(stat);
752  }
753 
754  /* --- [OPTIONAL] Node/Element number maps --- */
755  {
756  int *nmap = (int *)HECMW_malloc(sizeof(int) * n_node);
757  int *emap = (int *)HECMW_malloc(sizeof(int) * n_elem);
758  if (nmap == NULL || emap == NULL) {
759  HECMW_free(nmap); HECMW_free(emap);
760  HECMW_vis_print_exit("HECMW_malloc failed for node/elem maps");
761  }
762  /* Use global IDs if available, otherwise sequential */
763  if (mesh->global_node_ID != NULL) {
764  for (i = 0; i < n_node; i++) nmap[i] = mesh->global_node_ID[i];
765  } else {
766  for (i = 0; i < n_node; i++) nmap[i] = i + 1;
767  }
768  if (mesh->global_elem_ID != NULL) {
769  for (i = 0; i < n_elem; i++) emap[i] = mesh->global_elem_ID[i];
770  } else {
771  for (i = 0; i < n_elem; i++) emap[i] = i + 1;
772  }
773  stat = nc_put_var_int(ncid, var_node_num_map, nmap); NCERR(stat);
774  stat = nc_put_var_int(ncid, var_elem_num_map, emap); NCERR(stat);
775  HECMW_free(nmap); HECMW_free(emap);
776  }
777 
778  /* --- Nodal variable names --- */
779  if (total_nod_var > 0) {
780  char *name_buf = (char *)HECMW_malloc(total_nod_var * EXO_MAX_STR_LENGTH);
781  if (name_buf == NULL)
782  HECMW_vis_print_exit("HECMW_malloc failed for nodal variable names");
783  int vi = 0;
784  for (i = 0; i < data->nn_component; i++) {
785  const char **suffixes = get_comp_suffixes(data->nn_dof[i], data->node_label[i]);
786  for (k = 0; k < data->nn_dof[i]; k++) {
787  char vname[EXO_MAX_STR_LENGTH];
788  if (suffixes != NULL) {
789  snprintf(vname, sizeof(vname), "%s%s",
790  data->node_label[i], suffixes[k]);
791  } else {
792  snprintf(vname, sizeof(vname), "%s_%d",
793  data->node_label[i], k + 1);
794  }
795  pad_string(&name_buf[vi * EXO_MAX_STR_LENGTH], EXO_MAX_STR_LENGTH, vname);
796  vi++;
797  }
798  }
799  stat = nc_put_var_text(ncid, var_name_nod_var, name_buf);
800  NCERR(stat);
801  HECMW_free(name_buf);
802  }
803 
804  /* --- Element variable names --- */
805  if (total_elem_var > 0) {
806  char *name_buf = (char *)HECMW_malloc(total_elem_var * EXO_MAX_STR_LENGTH);
807  if (name_buf == NULL)
808  HECMW_vis_print_exit("HECMW_malloc failed for element variable names");
809  int vi = 0;
810  for (i = 0; i < data->ne_component; i++) {
811  const char **suffixes = get_comp_suffixes(data->ne_dof[i], data->elem_label[i]);
812  for (k = 0; k < data->ne_dof[i]; k++) {
813  char vname[EXO_MAX_STR_LENGTH];
814  if (suffixes != NULL) {
815  snprintf(vname, sizeof(vname), "%s%s",
816  data->elem_label[i], suffixes[k]);
817  } else {
818  snprintf(vname, sizeof(vname), "%s_%d",
819  data->elem_label[i], k + 1);
820  }
821  pad_string(&name_buf[vi * EXO_MAX_STR_LENGTH], EXO_MAX_STR_LENGTH, vname);
822  vi++;
823  }
824  }
825  stat = nc_put_var_text(ncid, var_name_elem_var, name_buf);
826  NCERR(stat);
827  HECMW_free(name_buf);
828  }
829 
830  /* --- Element variable truth table (all 1s: every var in every block) --- */
831  if (total_elem_var > 0) {
832  int *tab = (int *)HECMW_malloc(sizeof(int) * num_blocks * total_elem_var);
833  if (tab == NULL)
834  HECMW_vis_print_exit("HECMW_malloc failed for truth table");
835  for (i = 0; i < num_blocks * total_elem_var; i++) tab[i] = 1;
836  stat = nc_put_var_int(ncid, var_elem_var_tab, tab);
837  NCERR(stat);
838  HECMW_free(tab);
839  }
840 
841  } else {
842  /*====================================================================
843  * SUBSEQUENT CALLS (per_step==0 only): Open existing file for appending
844  *====================================================================*/
845  stat = nc_open(file_exo, NC_WRITE, &ncid);
846  NCERR(stat);
847  }
848 
849  /*========================================================================
850  * WRITE TIME-DEPENDENT DATA (every call)
851  *========================================================================*/
852 
853  /* --- Get time value from global data (look for TOTALTIME) --- */
854  double time_value = 0.0;
855  {
856  int gshift = 0;
857  for (i = 0; i < data->ng_component; i++) {
858  if (strcmp(data->global_label[i], "TOTALTIME") == 0) {
859  time_value = data->global_val_item[gshift];
860  break;
861  }
862  gshift += data->ng_dof[i];
863  }
864  }
865 
866  /* Write time value */
867  {
868  int varid;
869  /* per_step: each file has only 1 timestep (index=0)
870  * normal: timesteps accumulate in one file */
871  size_t ts_idx = per_step ? 0 : (size_t)time_step_count;
872  stat = nc_inq_varid(ncid, "time_whole", &varid); NCERR(stat);
873  stat = nc_put_var1_double(ncid, varid, &ts_idx, &time_value); NCERR(stat);
874  }
875 
876  /* --- Time step index for data arrays --- */
877  size_t ts_idx = per_step ? 0 : (size_t)time_step_count;
878 
879  /* --- Write nodal variables for this time step --- */
880  if (total_nod_var > 0) {
881  int vi = 0;
882  int comp_shift = 0;
883  for (i = 0; i < data->nn_component; i++) {
884  for (k = 0; k < data->nn_dof[i]; k++) {
885  char vname[EXO_MAX_STR_LENGTH];
886  snprintf(vname, sizeof(vname), "vals_nod_var%d", vi + 1);
887  int varid;
888  stat = nc_inq_varid(ncid, vname, &varid); NCERR(stat);
889 
890  double *vals = (double *)HECMW_malloc(sizeof(double) * n_node);
891  if (vals == NULL)
892  HECMW_vis_print_exit("HECMW_malloc failed for nodal values");
893  for (j = 0; j < n_node; j++) {
894  vals[j] = data->node_val_item[j * data_tot_n + comp_shift + k];
895  }
896  size_t start[2] = {ts_idx, 0};
897  size_t count[2] = {1, (size_t)n_node};
898  stat = nc_put_vara_double(ncid, varid, start, count, vals);
899  NCERR(stat);
900  HECMW_free(vals);
901  vi++;
902  }
903  comp_shift += data->nn_dof[i];
904  }
905  }
906 
907  /* --- Write element variables for this time step --- */
908  if (total_elem_var > 0) {
909  int vi = 0;
910  int comp_shift = 0;
911  for (i = 0; i < data->ne_component; i++) {
912  for (k = 0; k < data->ne_dof[i]; k++) {
913  /* Write for each block */
914  for (j = 0; j < num_blocks; j++) {
915  char vname[EXO_MAX_STR_LENGTH];
916  snprintf(vname, sizeof(vname), "vals_elem_var%deb%d", vi + 1, j + 1);
917  int varid;
918  stat = nc_inq_varid(ncid, vname, &varid); NCERR(stat);
919 
920  int ne = blocks[j].num_elem;
921  double *vals = (double *)HECMW_malloc(sizeof(double) * ne);
922  if (vals == NULL)
923  HECMW_vis_print_exit("HECMW_malloc failed for element values");
924  int ei_idx;
925  for (ei_idx = 0; ei_idx < ne; ei_idx++) {
926  int ei = blocks[j].elem_indices[ei_idx];
927  vals[ei_idx] = data->elem_val_item[ei * data_tot_e + comp_shift + k];
928  }
929  size_t start[2] = {ts_idx, 0};
930  size_t count[2] = {1, (size_t)ne};
931  stat = nc_put_vara_double(ncid, varid, start, count, vals);
932  NCERR(stat);
933  HECMW_free(vals);
934  }
935  vi++;
936  }
937  comp_shift += data->ne_dof[i];
938  }
939  }
940 
941  /*========================================================================
942  * CLEANUP
943  *========================================================================*/
944  stat = nc_close(ncid);
945  NCERR(stat);
946 
947  /* Free element block index arrays */
948  for (j = 0; j < num_blocks; j++) {
949  HECMW_free(blocks[j].elem_indices);
950  }
951 
952  time_step_count++;
953 }
954 
955 /*============================================================================
956  * Public interface (same pattern as HECMW_vtk_output / HECMW_bin_vtk_output)
957  *============================================================================*/
958 
959 /* output_type=18 (EXODUS): All timesteps in one file */
960 void HECMW_exodus_output(struct hecmwST_local_mesh *mesh,
961  struct hecmwST_result_data *data,
962  char *outfile, char *outfile1,
963  HECMW_Comm VIS_COMM)
964 {
965  exodus_output(mesh, data, outfile, outfile1, VIS_COMM, 0);
966 }
967 
968 /* output_type=19 (STEP_EXODUS): One file per timestep */
970  struct hecmwST_result_data *data,
971  char *outfile, char *outfile1,
972  HECMW_Comm VIS_COMM)
973 {
974  exodus_output(mesh, data, outfile, outfile1, VIS_COMM, 1);
975 }
976 
977 #else /* !WITH_NETCDF */
978 
979 #include <stdio.h>
980 #include "hecmw_struct.h"
981 #include "hecmw_result.h"
982 
983 /*============================================================================
984  * Stub functions when compiled without NetCDF support
985  *============================================================================*/
986 static void exodus_stub_warning(HECMW_Comm VIS_COMM)
987 {
988  int mynode;
989  HECMW_Comm_rank(VIS_COMM, &mynode);
990  if (mynode == 0) {
991  fprintf(stderr,
992  "WARNING: output_type=EXODUS requested but FrontISTR was not\n"
993  " compiled with NetCDF support.\n"
994  " Rebuild with -DWITH_NETCDF=ON (CMake) or\n"
995  " --with-netcdf (Makefile) to enable Exodus output.\n");
996  }
997 }
998 
1000  struct hecmwST_result_data *data,
1001  char *outfile, char *outfile1,
1002  HECMW_Comm VIS_COMM)
1003 {
1004  exodus_stub_warning(VIS_COMM);
1005  (void)mesh; (void)data; (void)outfile; (void)outfile1;
1006 }
1007 
1009  struct hecmwST_result_data *data,
1010  char *outfile, char *outfile1,
1011  HECMW_Comm VIS_COMM)
1012 {
1013  exodus_stub_warning(VIS_COMM);
1014  (void)mesh; (void)data; (void)outfile; (void)outfile1;
1015 }
1016 
1017 #endif /* WITH_NETCDF */
int HECMW_Comm_rank(HECMW_Comm comm, int *rank)
Definition: hecmw_comm.c:18
int HECMW_Comm_size(HECMW_Comm comm, int *size)
Definition: hecmw_comm.c:37
#define HECMW_FILENAME_LEN
Definition: hecmw_config.h:74
MPI_Comm HECMW_Comm
Definition: hecmw_config.h:30
int HECMW_ctrl_make_subdir(char *filename)
if(!(yy_init))
struct hecmwST_local_mesh * mesh
Definition: hecmw_repart.h:71
void HECMW_exodus_output(struct hecmwST_local_mesh *mesh, struct hecmwST_result_data *data, char *outfile, char *outfile1, HECMW_Comm VIS_COMM)
void HECMW_exodus_step_output(struct hecmwST_local_mesh *mesh, struct hecmwST_result_data *data, char *outfile, char *outfile1, HECMW_Comm VIS_COMM)
#define NULL
#define HECMW_free(ptr)
Definition: hecmw_malloc.h:24
#define HECMW_malloc(size)
Definition: hecmw_malloc.h:20
void HECMW_vis_print_exit(char *var)
integer(kind=kint) myrank
PARALLEL EXECUTION.
Definition: m_fstr.F90:99
void title()
Definition: neu2fstr.cpp:52
CNFData data
long long * elem_node_index
Definition: hecmw_struct.h:195