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 int 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  int jS = mesh->elem_node_index[elem_idx];
160  int jE = mesh->elem_node_index[elem_idx + 1];
161  int 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  int 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  strncpy(dir_exo, outfile1, HECMW_FILENAME_LEN - 1);
310  dir_exo[HECMW_FILENAME_LEN - 1] = '\0';
311  p = strrchr(dir_exo, '.');
312  if (p != NULL) *p = '\0'; /* "vis_out_psf.0001" -> "vis_out_psf" */
313  strncat(dir_exo, "_exo", HECMW_FILENAME_LEN - strlen(dir_exo) - 1);
314  /* dir_exo = "vis_out_psf_exo" */
315 
316  /* Build basename: strip directory prefix and timestep suffix from outfile */
317  strncpy(basename_exo, outfile, HECMW_FILENAME_LEN - 1);
318  basename_exo[HECMW_FILENAME_LEN - 1] = '\0';
319  p = strrchr(basename_exo, '/');
320  if (p != NULL) {
321  memmove(basename_exo, p + 1, strlen(p + 1) + 1);
322  }
323  p = strrchr(basename_exo, '.');
324  if (p != NULL) *p = '\0'; /* "vis_out_psf.0001" -> "vis_out_psf" */
325  }
326 
327  if (per_step) {
328  /*
329  * STEP_EXODUS mode: one file per timestep
330  * serial: {basename}.step{NNNN}.exo
331  * parallel: {basename}.step{NNNN}.exo.{nprocs}.{rank} (Nemesis)
332  */
333  if (petot == 1) {
334  snprintf(file_exo, HECMW_FILENAME_LEN, "%s/%s.step%04d.exo",
335  dir_exo, basename_exo, time_step_count);
336  } else {
337  snprintf(file_exo, HECMW_FILENAME_LEN, "%s/%s.step%04d.exo.%d.%d",
338  dir_exo, basename_exo, time_step_count, petot, myrank);
339  }
340  } else {
341  /*
342  * EXODUS mode: all timesteps in one file (only set path on first call)
343  * serial: {basename}.exo
344  * parallel: {basename}.exo.{nprocs}.{rank} (Nemesis)
345  */
346  if (time_step_count == 0) {
347  if (petot == 1) {
348  snprintf(file_exo, HECMW_FILENAME_LEN, "%s/%s.exo",
349  dir_exo, basename_exo);
350  } else {
351  snprintf(file_exo, HECMW_FILENAME_LEN, "%s/%s.exo.%d.%d",
352  dir_exo, basename_exo, petot, myrank);
353  }
354  }
355  /* For subsequent calls, file_exo is already set (static) */
356  }
357 
358  /* Create the directory (needed on first call, or every call for per_step) */
359  if (time_step_count == 0 || per_step) {
360  if (HECMW_ctrl_make_subdir(file_exo)) {
361  HECMW_vis_print_exit("ERROR: HEC-MW-VIS-E0009: Cannot open output directory for Exodus");
362  }
363  }
364 
365  /*========================================================================
366  * PASS 1: Scan elements and group into blocks by element type
367  *========================================================================*/
368  memset(blocks, 0, sizeof(blocks));
369  for (i = 0; i < n_elem; i++) {
370  int etype = mesh->elem_type[i];
371  int found = 0;
372  for (j = 0; j < num_blocks; j++) {
373  if (blocks[j].hecmw_type == etype) {
374  blocks[j].num_elem++;
375  found = 1;
376  break;
377  }
378  }
379  if (!found) {
380  if (num_blocks >= MAX_ELEM_BLOCKS) {
381  HECMW_vis_print_exit("ERROR: Too many element types for Exodus output");
382  }
383  const ElemTypeMap *map = get_elem_type_map(etype);
384  if (map == NULL) {
385  fprintf(stderr, "WARNING: Unknown element type %d, skipping\n", etype);
386  continue;
387  }
388  blocks[num_blocks].hecmw_type = etype;
389  blocks[num_blocks].exo_name = map->exo_name;
390  /* Determine nodes per element from first element of this type */
391  blocks[num_blocks].nod_per_elem = get_effective_n_nodes(mesh, i);
392  blocks[num_blocks].num_elem = 1;
393  num_blocks++;
394  }
395  }
396 
397  /* Allocate index arrays and fill element indices for each block */
398  for (j = 0; j < num_blocks; j++) {
399  blocks[j].elem_indices = (int *)HECMW_malloc(sizeof(int) * blocks[j].num_elem);
400  if (blocks[j].elem_indices == NULL)
401  HECMW_vis_print_exit("HECMW_malloc failed for block elem_indices");
402  blocks[j].num_elem = 0; /* reset count, will re-fill */
403  }
404  for (i = 0; i < n_elem; i++) {
405  int etype = mesh->elem_type[i];
406  for (j = 0; j < num_blocks; j++) {
407  if (blocks[j].hecmw_type == etype) {
408  blocks[j].elem_indices[blocks[j].num_elem] = i;
409  blocks[j].num_elem++;
410  break;
411  }
412  }
413  }
414 
415  /*========================================================================
416  * CREATE or OPEN file
417  *
418  * per_step==1: Always create a new file (one file per timestep)
419  * per_step==0: Create on first call, open for append on subsequent calls
420  *========================================================================*/
421  if (per_step || time_step_count == 0) {
422  /*====================================================================
423  * FIRST CALL: Create file and define full schema
424  *====================================================================*/
425  stat = nc_create(file_exo, NC_CLOBBER | NC_64BIT_OFFSET, &ncid);
426  NCERR(stat);
427  stat = nc_set_fill(ncid, NC_NOFILL, &old_fill_mode);
428  NCERR(stat);
429 
430  /*----------------------------------------------------------------
431  * Global Attributes [REQUIRED]
432  *----------------------------------------------------------------*/
433  {
434  float api_ver = EXO_API_VERSION;
435  float db_ver = EXO_DB_VERSION;
436  int word_sz = (int)sizeof(double);
437  int file_sz = EXO_FILE_SIZE;
438  char title[] = "FrontISTR Exodus output";
439 
440  stat = nc_put_att_float(ncid, NC_GLOBAL, "api_version", NC_FLOAT, 1, &api_ver);
441  NCERR(stat);
442  stat = nc_put_att_float(ncid, NC_GLOBAL, "version", NC_FLOAT, 1, &db_ver);
443  NCERR(stat);
444  stat = nc_put_att_int(ncid, NC_GLOBAL, "floating_point_word_size", NC_INT, 1, &word_sz);
445  NCERR(stat);
446  stat = nc_put_att_int(ncid, NC_GLOBAL, "file_size", NC_INT, 1, &file_sz);
447  NCERR(stat);
448  stat = nc_put_att_text(ncid, NC_GLOBAL, "title", strlen(title), title);
449  NCERR(stat);
450  }
451 
452  /*----------------------------------------------------------------
453  * Dimensions [REQUIRED]
454  *----------------------------------------------------------------*/
455  int dim_len_string, dim_len_name, dim_len_line, dim_four;
456  int dim_time_step, dim_num_dim, dim_num_nodes, dim_num_elem;
457  int dim_num_el_blk;
458  int dim_num_qa_rec;
459 
460  stat = nc_def_dim(ncid, "len_string", EXO_MAX_STR_LENGTH, &dim_len_string);
461  NCERR(stat);
462  stat = nc_def_dim(ncid, "len_name", EXO_MAX_STR_LENGTH, &dim_len_name);
463  NCERR(stat);
464  stat = nc_def_dim(ncid, "len_line", EXO_MAX_LINE_LENGTH, &dim_len_line);
465  NCERR(stat);
466  stat = nc_def_dim(ncid, "four", 4, &dim_four);
467  NCERR(stat);
468  stat = nc_def_dim(ncid, "time_step", NC_UNLIMITED, &dim_time_step);
469  NCERR(stat);
470  stat = nc_def_dim(ncid, "num_dim", 3, &dim_num_dim);
471  NCERR(stat);
472  stat = nc_def_dim(ncid, "num_nodes", (size_t)n_node, &dim_num_nodes);
473  NCERR(stat);
474  stat = nc_def_dim(ncid, "num_elem", (size_t)n_elem, &dim_num_elem);
475  NCERR(stat);
476  stat = nc_def_dim(ncid, "num_el_blk", (size_t)num_blocks, &dim_num_el_blk);
477  NCERR(stat);
478 
479  /* [OPTIONAL] QA record dimension */
480  stat = nc_def_dim(ncid, "num_qa_rec", 1, &dim_num_qa_rec);
481  NCERR(stat);
482 
483  /* Per-block dimensions */
484  int dim_el_in_blk[MAX_ELEM_BLOCKS];
485  int dim_nod_per_el[MAX_ELEM_BLOCKS];
486  for (j = 0; j < num_blocks; j++) {
487  char dname[EXO_MAX_STR_LENGTH];
488  snprintf(dname, EXO_MAX_STR_LENGTH, "num_el_in_blk%d", j + 1);
489  stat = nc_def_dim(ncid, dname, (size_t)blocks[j].num_elem, &dim_el_in_blk[j]);
490  NCERR(stat);
491  snprintf(dname, EXO_MAX_STR_LENGTH, "num_nod_per_el%d", j + 1);
492  stat = nc_def_dim(ncid, dname, (size_t)blocks[j].nod_per_elem, &dim_nod_per_el[j]);
493  NCERR(stat);
494  }
495 
496  /* Variable count dimensions */
497  int dim_num_nod_var = -1, dim_num_elem_var = -1;
498  if (total_nod_var > 0) {
499  stat = nc_def_dim(ncid, "num_nod_var", (size_t)total_nod_var, &dim_num_nod_var);
500  NCERR(stat);
501  }
502  if (total_elem_var > 0) {
503  stat = nc_def_dim(ncid, "num_elem_var", (size_t)total_elem_var, &dim_num_elem_var);
504  NCERR(stat);
505  }
506 
507  /*----------------------------------------------------------------
508  * Variables: Time [REQUIRED]
509  *----------------------------------------------------------------*/
510  int var_time_whole;
511  stat = nc_def_var(ncid, "time_whole", NC_DOUBLE, 1, &dim_time_step, &var_time_whole);
512  NCERR(stat);
513 
514  /*----------------------------------------------------------------
515  * Variables: Coordinates [REQUIRED]
516  *----------------------------------------------------------------*/
517  int var_coordx, var_coordy, var_coordz;
518  stat = nc_def_var(ncid, "coordx", NC_DOUBLE, 1, &dim_num_nodes, &var_coordx);
519  NCERR(stat);
520  stat = nc_def_var(ncid, "coordy", NC_DOUBLE, 1, &dim_num_nodes, &var_coordy);
521  NCERR(stat);
522  stat = nc_def_var(ncid, "coordz", NC_DOUBLE, 1, &dim_num_nodes, &var_coordz);
523  NCERR(stat);
524 
525  /* coor_names [REQUIRED] */
526  int var_coor_names;
527  {
528  int dims2[2] = {dim_num_dim, dim_len_name};
529  stat = nc_def_var(ncid, "coor_names", NC_CHAR, 2, dims2, &var_coor_names);
530  NCERR(stat);
531  }
532 
533  /*----------------------------------------------------------------
534  * Variables: Element Block info [REQUIRED]
535  *----------------------------------------------------------------*/
536  int var_eb_status, var_eb_prop1;
537  stat = nc_def_var(ncid, "eb_status", NC_INT, 1, &dim_num_el_blk, &var_eb_status);
538  NCERR(stat);
539  stat = nc_def_var(ncid, "eb_prop1", NC_INT, 1, &dim_num_el_blk, &var_eb_prop1);
540  NCERR(stat);
541  stat = nc_put_att_text(ncid, var_eb_prop1, "name", 2, "ID");
542  NCERR(stat);
543 
544  /* [OPTIONAL] eb_names */
545  int var_eb_names;
546  {
547  int dims2[2] = {dim_num_el_blk, dim_len_name};
548  stat = nc_def_var(ncid, "eb_names", NC_CHAR, 2, dims2, &var_eb_names);
549  NCERR(stat);
550  }
551 
552  /* Connectivity per block [REQUIRED] */
553  int var_connect[MAX_ELEM_BLOCKS];
554  for (j = 0; j < num_blocks; j++) {
555  char vname[EXO_MAX_STR_LENGTH];
556  snprintf(vname, EXO_MAX_STR_LENGTH, "connect%d", j + 1);
557  int dims2[2] = {dim_el_in_blk[j], dim_nod_per_el[j]};
558  stat = nc_def_var(ncid, vname, NC_INT, 2, dims2, &var_connect[j]);
559  NCERR(stat);
560  stat = nc_put_att_text(ncid, var_connect[j], "elem_type",
561  strlen(blocks[j].exo_name), blocks[j].exo_name);
562  NCERR(stat);
563  }
564 
565  /*----------------------------------------------------------------
566  * Variables: QA records [OPTIONAL - nice to have]
567  *----------------------------------------------------------------*/
568  int var_qa_records;
569  {
570  int dims3[3] = {dim_num_qa_rec, dim_four, dim_len_string};
571  stat = nc_def_var(ncid, "qa_records", NC_CHAR, 3, dims3, &var_qa_records);
572  NCERR(stat);
573  }
574 
575  /*----------------------------------------------------------------
576  * Variables: ID maps [OPTIONAL - nice to have for parallel]
577  *----------------------------------------------------------------*/
578  int var_node_num_map, var_elem_num_map;
579  stat = nc_def_var(ncid, "node_num_map", NC_INT, 1, &dim_num_nodes, &var_node_num_map);
580  NCERR(stat);
581  stat = nc_def_var(ncid, "elem_num_map", NC_INT, 1, &dim_num_elem, &var_elem_num_map);
582  NCERR(stat);
583 
584  /*----------------------------------------------------------------
585  * Variables: Nodal result variable names and data [REQUIRED]
586  *----------------------------------------------------------------*/
587  int var_name_nod_var = -1;
588  int var_vals_nod[512]; /* enough for many variables */
589  if (total_nod_var > 0) {
590  int dims2[2] = {dim_num_nod_var, dim_len_name};
591  stat = nc_def_var(ncid, "name_nod_var", NC_CHAR, 2, dims2, &var_name_nod_var);
592  NCERR(stat);
593 
594  for (i = 0; i < total_nod_var; i++) {
595  char vname[EXO_MAX_STR_LENGTH];
596  snprintf(vname, EXO_MAX_STR_LENGTH, "vals_nod_var%d", i + 1);
597  int dims_ts_nn[2] = {dim_time_step, dim_num_nodes};
598  stat = nc_def_var(ncid, vname, NC_DOUBLE, 2, dims_ts_nn, &var_vals_nod[i]);
599  NCERR(stat);
600  }
601  }
602 
603  /*----------------------------------------------------------------
604  * Variables: Element result variable names and data [REQUIRED]
605  *----------------------------------------------------------------*/
606  int var_name_elem_var = -1;
607  int var_vals_elem[512]; /* [var_idx * num_blocks + blk_idx] */
608  int var_elem_var_tab = -1;
609  if (total_elem_var > 0) {
610  int dims2[2] = {dim_num_elem_var, dim_len_name};
611  stat = nc_def_var(ncid, "name_elem_var", NC_CHAR, 2, dims2, &var_name_elem_var);
612  NCERR(stat);
613 
614  /* Truth table */
615  {
616  int dims_tt[2] = {dim_num_el_blk, dim_num_elem_var};
617  stat = nc_def_var(ncid, "elem_var_tab", NC_INT, 2, dims_tt, &var_elem_var_tab);
618  NCERR(stat);
619  }
620 
621  /* vals_elem_var{v}eb{b} for each variable and block */
622  for (i = 0; i < total_elem_var; i++) {
623  for (j = 0; j < num_blocks; j++) {
624  char vname[EXO_MAX_STR_LENGTH];
625  snprintf(vname, EXO_MAX_STR_LENGTH, "vals_elem_var%deb%d", i + 1, j + 1);
626  int dims_ts_eb[2] = {dim_time_step, dim_el_in_blk[j]};
627  stat = nc_def_var(ncid, vname, NC_DOUBLE, 2, dims_ts_eb,
628  &var_vals_elem[i * num_blocks + j]);
629  NCERR(stat);
630  }
631  }
632  }
633 
634  /*================================================================
635  * END DEFINE MODE
636  *================================================================*/
637  stat = nc_enddef(ncid);
638  NCERR(stat);
639 
640  /*================================================================
641  * WRITE STATIC DATA (coordinates, connectivity, maps, etc.)
642  * These are written only once.
643  *================================================================*/
644 
645  /* --- Coordinates (XYZ separated) --- */
646  {
647  double *cx = (double *)HECMW_malloc(sizeof(double) * n_node);
648  double *cy = (double *)HECMW_malloc(sizeof(double) * n_node);
649  double *cz = (double *)HECMW_malloc(sizeof(double) * n_node);
650  if (cx == NULL || cy == NULL || cz == NULL) {
651  HECMW_free(cx); HECMW_free(cy); HECMW_free(cz);
652  HECMW_vis_print_exit("HECMW_malloc failed for coordinate arrays");
653  }
654  for (i = 0; i < n_node; i++) {
655  cx[i] = mesh->node[3 * i];
656  cy[i] = mesh->node[3 * i + 1];
657  cz[i] = mesh->node[3 * i + 2];
658  }
659  stat = nc_put_var_double(ncid, var_coordx, cx); NCERR(stat);
660  stat = nc_put_var_double(ncid, var_coordy, cy); NCERR(stat);
661  stat = nc_put_var_double(ncid, var_coordz, cz); NCERR(stat);
662  HECMW_free(cx); HECMW_free(cy); HECMW_free(cz);
663  }
664 
665  /* --- Coordinate names --- */
666  {
667  char names[3][EXO_MAX_STR_LENGTH];
668  pad_string(names[0], EXO_MAX_STR_LENGTH, "x");
669  pad_string(names[1], EXO_MAX_STR_LENGTH, "y");
670  pad_string(names[2], EXO_MAX_STR_LENGTH, "z");
671  stat = nc_put_var_text(ncid, var_coor_names, &names[0][0]);
672  NCERR(stat);
673  }
674 
675  /* --- Element block status and IDs --- */
676  {
677  int *eb_stat = (int *)HECMW_malloc(sizeof(int) * num_blocks);
678  int *eb_ids = (int *)HECMW_malloc(sizeof(int) * num_blocks);
679  if (eb_stat == NULL || eb_ids == NULL) {
680  HECMW_free(eb_stat); HECMW_free(eb_ids);
681  HECMW_vis_print_exit("HECMW_malloc failed for element block arrays");
682  }
683  for (j = 0; j < num_blocks; j++) {
684  eb_stat[j] = 1; /* active */
685  eb_ids[j] = j + 1; /* block ID = 1,2,3,... */
686  }
687  stat = nc_put_var_int(ncid, var_eb_status, eb_stat); NCERR(stat);
688  stat = nc_put_var_int(ncid, var_eb_prop1, eb_ids); NCERR(stat);
689  HECMW_free(eb_stat); HECMW_free(eb_ids);
690  }
691 
692  /* --- [OPTIONAL] Element block names --- */
693  {
694  char *names_buf = (char *)HECMW_malloc(num_blocks * EXO_MAX_STR_LENGTH);
695  if (names_buf == NULL)
696  HECMW_vis_print_exit("HECMW_malloc failed for block names");
697  for (j = 0; j < num_blocks; j++) {
698  char label[EXO_MAX_STR_LENGTH];
699  snprintf(label, EXO_MAX_STR_LENGTH, "block_%d_%s", j + 1, blocks[j].exo_name);
700  pad_string(&names_buf[j * EXO_MAX_STR_LENGTH], EXO_MAX_STR_LENGTH, label);
701  }
702  stat = nc_put_var_text(ncid, var_eb_names, names_buf);
703  NCERR(stat);
704  HECMW_free(names_buf);
705  }
706 
707  /* --- Connectivity per block --- */
708  for (j = 0; j < num_blocks; j++) {
709  int ne = blocks[j].num_elem;
710  int nn = blocks[j].nod_per_elem;
711  int etype = blocks[j].hecmw_type;
712  int node_shift = get_node_shift(etype);
713  int *conn = (int *)HECMW_malloc(sizeof(int) * ne * nn);
714  if (conn == NULL)
715  HECMW_vis_print_exit("HECMW_malloc failed for connectivity");
716 
717  for (i = 0; i < ne; i++) {
718  int ei = blocks[j].elem_indices[i]; /* original element index */
719  jS = mesh->elem_node_index[ei];
720  jE = mesh->elem_node_index[ei + 1];
721 
722  if (etype == 342) {
723  /* TET10: reorder nodes (same table as VTK output) */
724  for (k = 0; k < nn; k++) {
725  conn[i * nn + k] = mesh->elem_node_item[jS + table342[k]];
726  /* NOTE: Exodus uses 1-based indexing (no -1) */
727  }
728  } else {
729  for (k = 0; k < nn; k++) {
730  conn[i * nn + k] = mesh->elem_node_item[jS + k];
731  }
732  }
733  }
734  stat = nc_put_var_int(ncid, var_connect[j], conn);
735  NCERR(stat);
736  HECMW_free(conn);
737  }
738 
739  /* --- [OPTIONAL] QA records --- */
740  {
741  char qa_data[1][4][EXO_MAX_STR_LENGTH];
742  time_t now = time(NULL);
743  struct tm *t = localtime(&now);
744  char datestr[32], timestr[32];
745  strftime(datestr, sizeof(datestr), "%Y%m%d", t);
746  strftime(timestr, sizeof(timestr), "%H:%M:%S", t);
747 
748  pad_string(qa_data[0][0], EXO_MAX_STR_LENGTH, "FrontISTR");
749  pad_string(qa_data[0][1], EXO_MAX_STR_LENGTH, "1.0");
750  pad_string(qa_data[0][2], EXO_MAX_STR_LENGTH, datestr);
751  pad_string(qa_data[0][3], EXO_MAX_STR_LENGTH, timestr);
752  stat = nc_put_var_text(ncid, var_qa_records, &qa_data[0][0][0]);
753  NCERR(stat);
754  }
755 
756  /* --- [OPTIONAL] Node/Element number maps --- */
757  {
758  int *nmap = (int *)HECMW_malloc(sizeof(int) * n_node);
759  int *emap = (int *)HECMW_malloc(sizeof(int) * n_elem);
760  if (nmap == NULL || emap == NULL) {
761  HECMW_free(nmap); HECMW_free(emap);
762  HECMW_vis_print_exit("HECMW_malloc failed for node/elem maps");
763  }
764  /* Use global IDs if available, otherwise sequential */
765  if (mesh->global_node_ID != NULL) {
766  for (i = 0; i < n_node; i++) nmap[i] = mesh->global_node_ID[i];
767  } else {
768  for (i = 0; i < n_node; i++) nmap[i] = i + 1;
769  }
770  if (mesh->global_elem_ID != NULL) {
771  for (i = 0; i < n_elem; i++) emap[i] = mesh->global_elem_ID[i];
772  } else {
773  for (i = 0; i < n_elem; i++) emap[i] = i + 1;
774  }
775  stat = nc_put_var_int(ncid, var_node_num_map, nmap); NCERR(stat);
776  stat = nc_put_var_int(ncid, var_elem_num_map, emap); NCERR(stat);
777  HECMW_free(nmap); HECMW_free(emap);
778  }
779 
780  /* --- Nodal variable names --- */
781  if (total_nod_var > 0) {
782  char *name_buf = (char *)HECMW_malloc(total_nod_var * EXO_MAX_STR_LENGTH);
783  if (name_buf == NULL)
784  HECMW_vis_print_exit("HECMW_malloc failed for nodal variable names");
785  int vi = 0;
786  for (i = 0; i < data->nn_component; i++) {
787  const char **suffixes = get_comp_suffixes(data->nn_dof[i], data->node_label[i]);
788  for (k = 0; k < data->nn_dof[i]; k++) {
789  char vname[EXO_MAX_STR_LENGTH];
790  if (suffixes != NULL) {
791  snprintf(vname, EXO_MAX_STR_LENGTH, "%s%s",
792  data->node_label[i], suffixes[k]);
793  } else {
794  snprintf(vname, EXO_MAX_STR_LENGTH, "%s_%d",
795  data->node_label[i], k + 1);
796  }
797  pad_string(&name_buf[vi * EXO_MAX_STR_LENGTH], EXO_MAX_STR_LENGTH, vname);
798  vi++;
799  }
800  }
801  stat = nc_put_var_text(ncid, var_name_nod_var, name_buf);
802  NCERR(stat);
803  HECMW_free(name_buf);
804  }
805 
806  /* --- Element variable names --- */
807  if (total_elem_var > 0) {
808  char *name_buf = (char *)HECMW_malloc(total_elem_var * EXO_MAX_STR_LENGTH);
809  if (name_buf == NULL)
810  HECMW_vis_print_exit("HECMW_malloc failed for element variable names");
811  int vi = 0;
812  for (i = 0; i < data->ne_component; i++) {
813  const char **suffixes = get_comp_suffixes(data->ne_dof[i], data->elem_label[i]);
814  for (k = 0; k < data->ne_dof[i]; k++) {
815  char vname[EXO_MAX_STR_LENGTH];
816  if (suffixes != NULL) {
817  snprintf(vname, EXO_MAX_STR_LENGTH, "%s%s",
818  data->elem_label[i], suffixes[k]);
819  } else {
820  snprintf(vname, EXO_MAX_STR_LENGTH, "%s_%d",
821  data->elem_label[i], k + 1);
822  }
823  pad_string(&name_buf[vi * EXO_MAX_STR_LENGTH], EXO_MAX_STR_LENGTH, vname);
824  vi++;
825  }
826  }
827  stat = nc_put_var_text(ncid, var_name_elem_var, name_buf);
828  NCERR(stat);
829  HECMW_free(name_buf);
830  }
831 
832  /* --- Element variable truth table (all 1s: every var in every block) --- */
833  if (total_elem_var > 0) {
834  int *tab = (int *)HECMW_malloc(sizeof(int) * num_blocks * total_elem_var);
835  if (tab == NULL)
836  HECMW_vis_print_exit("HECMW_malloc failed for truth table");
837  for (i = 0; i < num_blocks * total_elem_var; i++) tab[i] = 1;
838  stat = nc_put_var_int(ncid, var_elem_var_tab, tab);
839  NCERR(stat);
840  HECMW_free(tab);
841  }
842 
843  } else {
844  /*====================================================================
845  * SUBSEQUENT CALLS (per_step==0 only): Open existing file for appending
846  *====================================================================*/
847  stat = nc_open(file_exo, NC_WRITE, &ncid);
848  NCERR(stat);
849  }
850 
851  /*========================================================================
852  * WRITE TIME-DEPENDENT DATA (every call)
853  *========================================================================*/
854 
855  /* --- Get time value from global data (look for TOTALTIME) --- */
856  double time_value = 0.0;
857  {
858  int gshift = 0;
859  for (i = 0; i < data->ng_component; i++) {
860  if (strcmp(data->global_label[i], "TOTALTIME") == 0) {
861  time_value = data->global_val_item[gshift];
862  break;
863  }
864  gshift += data->ng_dof[i];
865  }
866  }
867 
868  /* Write time value */
869  {
870  int varid;
871  /* per_step: each file has only 1 timestep (index=0)
872  * normal: timesteps accumulate in one file */
873  size_t ts_idx = per_step ? 0 : (size_t)time_step_count;
874  stat = nc_inq_varid(ncid, "time_whole", &varid); NCERR(stat);
875  stat = nc_put_var1_double(ncid, varid, &ts_idx, &time_value); NCERR(stat);
876  }
877 
878  /* --- Time step index for data arrays --- */
879  size_t ts_idx = per_step ? 0 : (size_t)time_step_count;
880 
881  /* --- Write nodal variables for this time step --- */
882  if (total_nod_var > 0) {
883  int vi = 0;
884  int comp_shift = 0;
885  for (i = 0; i < data->nn_component; i++) {
886  for (k = 0; k < data->nn_dof[i]; k++) {
887  char vname[EXO_MAX_STR_LENGTH];
888  snprintf(vname, EXO_MAX_STR_LENGTH, "vals_nod_var%d", vi + 1);
889  int varid;
890  stat = nc_inq_varid(ncid, vname, &varid); NCERR(stat);
891 
892  double *vals = (double *)HECMW_malloc(sizeof(double) * n_node);
893  if (vals == NULL)
894  HECMW_vis_print_exit("HECMW_malloc failed for nodal values");
895  for (j = 0; j < n_node; j++) {
896  vals[j] = data->node_val_item[j * data_tot_n + comp_shift + k];
897  }
898  size_t start[2] = {ts_idx, 0};
899  size_t count[2] = {1, (size_t)n_node};
900  stat = nc_put_vara_double(ncid, varid, start, count, vals);
901  NCERR(stat);
902  HECMW_free(vals);
903  vi++;
904  }
905  comp_shift += data->nn_dof[i];
906  }
907  }
908 
909  /* --- Write element variables for this time step --- */
910  if (total_elem_var > 0) {
911  int vi = 0;
912  int comp_shift = 0;
913  for (i = 0; i < data->ne_component; i++) {
914  for (k = 0; k < data->ne_dof[i]; k++) {
915  /* Write for each block */
916  for (j = 0; j < num_blocks; j++) {
917  char vname[EXO_MAX_STR_LENGTH];
918  snprintf(vname, EXO_MAX_STR_LENGTH, "vals_elem_var%deb%d", vi + 1, j + 1);
919  int varid;
920  stat = nc_inq_varid(ncid, vname, &varid); NCERR(stat);
921 
922  int ne = blocks[j].num_elem;
923  double *vals = (double *)HECMW_malloc(sizeof(double) * ne);
924  if (vals == NULL)
925  HECMW_vis_print_exit("HECMW_malloc failed for element values");
926  int ei_idx;
927  for (ei_idx = 0; ei_idx < ne; ei_idx++) {
928  int ei = blocks[j].elem_indices[ei_idx];
929  vals[ei_idx] = data->elem_val_item[ei * data_tot_e + comp_shift + k];
930  }
931  size_t start[2] = {ts_idx, 0};
932  size_t count[2] = {1, (size_t)ne};
933  stat = nc_put_vara_double(ncid, varid, start, count, vals);
934  NCERR(stat);
935  HECMW_free(vals);
936  }
937  vi++;
938  }
939  comp_shift += data->ne_dof[i];
940  }
941  }
942 
943  /*========================================================================
944  * CLEANUP
945  *========================================================================*/
946  stat = nc_close(ncid);
947  NCERR(stat);
948 
949  /* Free element block index arrays */
950  for (j = 0; j < num_blocks; j++) {
951  HECMW_free(blocks[j].elem_indices);
952  }
953 
954  time_step_count++;
955 }
956 
957 /*============================================================================
958  * Public interface (same pattern as HECMW_vtk_output / HECMW_bin_vtk_output)
959  *============================================================================*/
960 
961 /* output_type=18 (EXODUS): All timesteps in one file */
962 void HECMW_exodus_output(struct hecmwST_local_mesh *mesh,
963  struct hecmwST_result_data *data,
964  char *outfile, char *outfile1,
965  HECMW_Comm VIS_COMM)
966 {
967  exodus_output(mesh, data, outfile, outfile1, VIS_COMM, 0);
968 }
969 
970 /* output_type=19 (STEP_EXODUS): One file per timestep */
972  struct hecmwST_result_data *data,
973  char *outfile, char *outfile1,
974  HECMW_Comm VIS_COMM)
975 {
976  exodus_output(mesh, data, outfile, outfile1, VIS_COMM, 1);
977 }
978 
979 #else /* !WITH_NETCDF */
980 
981 #include <stdio.h>
982 #include "hecmw_struct.h"
983 #include "hecmw_result.h"
984 
985 /*============================================================================
986  * Stub functions when compiled without NetCDF support
987  *============================================================================*/
988 static void exodus_stub_warning(HECMW_Comm VIS_COMM)
989 {
990  int mynode;
991  HECMW_Comm_rank(VIS_COMM, &mynode);
992  if (mynode == 0) {
993  fprintf(stderr,
994  "WARNING: output_type=EXODUS requested but FrontISTR was not\n"
995  " compiled with NetCDF support.\n"
996  " Rebuild with -DWITH_NETCDF=ON (CMake) or\n"
997  " --with-netcdf (Makefile) to enable Exodus output.\n");
998  }
999 }
1000 
1002  struct hecmwST_result_data *data,
1003  char *outfile, char *outfile1,
1004  HECMW_Comm VIS_COMM)
1005 {
1006  exodus_stub_warning(VIS_COMM);
1007  (void)mesh; (void)data; (void)outfile; (void)outfile1;
1008 }
1009 
1011  struct hecmwST_result_data *data,
1012  char *outfile, char *outfile1,
1013  HECMW_Comm VIS_COMM)
1014 {
1015  exodus_stub_warning(VIS_COMM);
1016  (void)mesh; (void)data; (void)outfile; (void)outfile1;
1017 }
1018 
1019 #endif /* WITH_NETCDF */
if(!(yy_init))
Definition: hecmw_ablex.c:1823
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)
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:97
void title()
Definition: neu2fstr.cpp:52
CNFData data
long long * elem_node_index
Definition: hecmw_struct.h:195