FrontISTR  5.8.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  *****************************************************************************/
7 #include "hecmw_fstr_output_exo.h"
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  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, EXO_MAX_STR_LENGTH, "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, EXO_MAX_STR_LENGTH, "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, EXO_MAX_STR_LENGTH, "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, EXO_MAX_STR_LENGTH, "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, EXO_MAX_STR_LENGTH, "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  for (i = 0; i < n_node; i++) {
649  cx[i] = mesh->node[3 * i];
650  cy[i] = mesh->node[3 * i + 1];
651  cz[i] = mesh->node[3 * i + 2];
652  }
653  stat = nc_put_var_double(ncid, var_coordx, cx); NCERR(stat);
654  stat = nc_put_var_double(ncid, var_coordy, cy); NCERR(stat);
655  stat = nc_put_var_double(ncid, var_coordz, cz); NCERR(stat);
656  HECMW_free(cx); HECMW_free(cy); HECMW_free(cz);
657  }
658 
659  /* --- Coordinate names --- */
660  {
661  char names[3][EXO_MAX_STR_LENGTH];
662  pad_string(names[0], EXO_MAX_STR_LENGTH, "x");
663  pad_string(names[1], EXO_MAX_STR_LENGTH, "y");
664  pad_string(names[2], EXO_MAX_STR_LENGTH, "z");
665  stat = nc_put_var_text(ncid, var_coor_names, &names[0][0]);
666  NCERR(stat);
667  }
668 
669  /* --- Element block status and IDs --- */
670  {
671  int *eb_stat = (int *)HECMW_malloc(sizeof(int) * num_blocks);
672  int *eb_ids = (int *)HECMW_malloc(sizeof(int) * num_blocks);
673  for (j = 0; j < num_blocks; j++) {
674  eb_stat[j] = 1; /* active */
675  eb_ids[j] = j + 1; /* block ID = 1,2,3,... */
676  }
677  stat = nc_put_var_int(ncid, var_eb_status, eb_stat); NCERR(stat);
678  stat = nc_put_var_int(ncid, var_eb_prop1, eb_ids); NCERR(stat);
679  HECMW_free(eb_stat); HECMW_free(eb_ids);
680  }
681 
682  /* --- [OPTIONAL] Element block names --- */
683  {
684  char *names_buf = (char *)HECMW_malloc(num_blocks * EXO_MAX_STR_LENGTH);
685  for (j = 0; j < num_blocks; j++) {
686  char label[EXO_MAX_STR_LENGTH];
687  snprintf(label, EXO_MAX_STR_LENGTH, "block_%d_%s", j + 1, blocks[j].exo_name);
688  pad_string(&names_buf[j * EXO_MAX_STR_LENGTH], EXO_MAX_STR_LENGTH, label);
689  }
690  stat = nc_put_var_text(ncid, var_eb_names, names_buf);
691  NCERR(stat);
692  HECMW_free(names_buf);
693  }
694 
695  /* --- Connectivity per block --- */
696  for (j = 0; j < num_blocks; j++) {
697  int ne = blocks[j].num_elem;
698  int nn = blocks[j].nod_per_elem;
699  int etype = blocks[j].hecmw_type;
700  int node_shift = get_node_shift(etype);
701  int *conn = (int *)HECMW_malloc(sizeof(int) * ne * nn);
702 
703  for (i = 0; i < ne; i++) {
704  int ei = blocks[j].elem_indices[i]; /* original element index */
705  jS = mesh->elem_node_index[ei];
706  jE = mesh->elem_node_index[ei + 1];
707 
708  if (etype == 342) {
709  /* TET10: reorder nodes (same table as VTK output) */
710  for (k = 0; k < nn; k++) {
711  conn[i * nn + k] = mesh->elem_node_item[jS + table342[k]];
712  /* NOTE: Exodus uses 1-based indexing (no -1) */
713  }
714  } else {
715  for (k = 0; k < nn; k++) {
716  conn[i * nn + k] = mesh->elem_node_item[jS + k];
717  }
718  }
719  }
720  stat = nc_put_var_int(ncid, var_connect[j], conn);
721  NCERR(stat);
722  HECMW_free(conn);
723  }
724 
725  /* --- [OPTIONAL] QA records --- */
726  {
727  char qa_data[1][4][EXO_MAX_STR_LENGTH];
728  time_t now = time(NULL);
729  struct tm *t = localtime(&now);
730  char datestr[32], timestr[32];
731  strftime(datestr, sizeof(datestr), "%Y%m%d", t);
732  strftime(timestr, sizeof(timestr), "%H:%M:%S", t);
733 
734  pad_string(qa_data[0][0], EXO_MAX_STR_LENGTH, "FrontISTR");
735  pad_string(qa_data[0][1], EXO_MAX_STR_LENGTH, "1.0");
736  pad_string(qa_data[0][2], EXO_MAX_STR_LENGTH, datestr);
737  pad_string(qa_data[0][3], EXO_MAX_STR_LENGTH, timestr);
738  stat = nc_put_var_text(ncid, var_qa_records, &qa_data[0][0][0]);
739  NCERR(stat);
740  }
741 
742  /* --- [OPTIONAL] Node/Element number maps --- */
743  {
744  int *nmap = (int *)HECMW_malloc(sizeof(int) * n_node);
745  int *emap = (int *)HECMW_malloc(sizeof(int) * n_elem);
746  /* Use global IDs if available, otherwise sequential */
747  if (mesh->global_node_ID != NULL) {
748  for (i = 0; i < n_node; i++) nmap[i] = mesh->global_node_ID[i];
749  } else {
750  for (i = 0; i < n_node; i++) nmap[i] = i + 1;
751  }
752  if (mesh->global_elem_ID != NULL) {
753  for (i = 0; i < n_elem; i++) emap[i] = mesh->global_elem_ID[i];
754  } else {
755  for (i = 0; i < n_elem; i++) emap[i] = i + 1;
756  }
757  stat = nc_put_var_int(ncid, var_node_num_map, nmap); NCERR(stat);
758  stat = nc_put_var_int(ncid, var_elem_num_map, emap); NCERR(stat);
759  HECMW_free(nmap); HECMW_free(emap);
760  }
761 
762  /* --- Nodal variable names --- */
763  if (total_nod_var > 0) {
764  char *name_buf = (char *)HECMW_malloc(total_nod_var * EXO_MAX_STR_LENGTH);
765  int vi = 0;
766  for (i = 0; i < data->nn_component; i++) {
767  const char **suffixes = get_comp_suffixes(data->nn_dof[i], data->node_label[i]);
768  for (k = 0; k < data->nn_dof[i]; k++) {
769  char vname[EXO_MAX_STR_LENGTH];
770  if (suffixes != NULL) {
771  snprintf(vname, EXO_MAX_STR_LENGTH, "%s%s",
772  data->node_label[i], suffixes[k]);
773  } else {
774  snprintf(vname, EXO_MAX_STR_LENGTH, "%s_%d",
775  data->node_label[i], k + 1);
776  }
777  pad_string(&name_buf[vi * EXO_MAX_STR_LENGTH], EXO_MAX_STR_LENGTH, vname);
778  vi++;
779  }
780  }
781  stat = nc_put_var_text(ncid, var_name_nod_var, name_buf);
782  NCERR(stat);
783  HECMW_free(name_buf);
784  }
785 
786  /* --- Element variable names --- */
787  if (total_elem_var > 0) {
788  char *name_buf = (char *)HECMW_malloc(total_elem_var * EXO_MAX_STR_LENGTH);
789  int vi = 0;
790  for (i = 0; i < data->ne_component; i++) {
791  const char **suffixes = get_comp_suffixes(data->ne_dof[i], data->elem_label[i]);
792  for (k = 0; k < data->ne_dof[i]; k++) {
793  char vname[EXO_MAX_STR_LENGTH];
794  if (suffixes != NULL) {
795  snprintf(vname, EXO_MAX_STR_LENGTH, "%s%s",
796  data->elem_label[i], suffixes[k]);
797  } else {
798  snprintf(vname, EXO_MAX_STR_LENGTH, "%s_%d",
799  data->elem_label[i], k + 1);
800  }
801  pad_string(&name_buf[vi * EXO_MAX_STR_LENGTH], EXO_MAX_STR_LENGTH, vname);
802  vi++;
803  }
804  }
805  stat = nc_put_var_text(ncid, var_name_elem_var, name_buf);
806  NCERR(stat);
807  HECMW_free(name_buf);
808  }
809 
810  /* --- Element variable truth table (all 1s: every var in every block) --- */
811  if (total_elem_var > 0) {
812  int *tab = (int *)HECMW_malloc(sizeof(int) * num_blocks * total_elem_var);
813  for (i = 0; i < num_blocks * total_elem_var; i++) tab[i] = 1;
814  stat = nc_put_var_int(ncid, var_elem_var_tab, tab);
815  NCERR(stat);
816  HECMW_free(tab);
817  }
818 
819  } else {
820  /*====================================================================
821  * SUBSEQUENT CALLS (per_step==0 only): Open existing file for appending
822  *====================================================================*/
823  stat = nc_open(file_exo, NC_WRITE, &ncid);
824  NCERR(stat);
825  }
826 
827  /*========================================================================
828  * WRITE TIME-DEPENDENT DATA (every call)
829  *========================================================================*/
830 
831  /* --- Get time value from global data (look for TOTALTIME) --- */
832  double time_value = 0.0;
833  {
834  int gshift = 0;
835  for (i = 0; i < data->ng_component; i++) {
836  if (strcmp(data->global_label[i], "TOTALTIME") == 0) {
837  time_value = data->global_val_item[gshift];
838  break;
839  }
840  gshift += data->ng_dof[i];
841  }
842  }
843 
844  /* Write time value */
845  {
846  int varid;
847  /* per_step: each file has only 1 timestep (index=0)
848  * normal: timesteps accumulate in one file */
849  size_t ts_idx = per_step ? 0 : (size_t)time_step_count;
850  stat = nc_inq_varid(ncid, "time_whole", &varid); NCERR(stat);
851  stat = nc_put_var1_double(ncid, varid, &ts_idx, &time_value); NCERR(stat);
852  }
853 
854  /* --- Time step index for data arrays --- */
855  size_t ts_idx = per_step ? 0 : (size_t)time_step_count;
856 
857  /* --- Write nodal variables for this time step --- */
858  if (total_nod_var > 0) {
859  int vi = 0;
860  int comp_shift = 0;
861  for (i = 0; i < data->nn_component; i++) {
862  for (k = 0; k < data->nn_dof[i]; k++) {
863  char vname[EXO_MAX_STR_LENGTH];
864  snprintf(vname, EXO_MAX_STR_LENGTH, "vals_nod_var%d", vi + 1);
865  int varid;
866  stat = nc_inq_varid(ncid, vname, &varid); NCERR(stat);
867 
868  double *vals = (double *)HECMW_malloc(sizeof(double) * n_node);
869  for (j = 0; j < n_node; j++) {
870  vals[j] = data->node_val_item[j * data_tot_n + comp_shift + k];
871  }
872  size_t start[2] = {ts_idx, 0};
873  size_t count[2] = {1, (size_t)n_node};
874  stat = nc_put_vara_double(ncid, varid, start, count, vals);
875  NCERR(stat);
876  HECMW_free(vals);
877  vi++;
878  }
879  comp_shift += data->nn_dof[i];
880  }
881  }
882 
883  /* --- Write element variables for this time step --- */
884  if (total_elem_var > 0) {
885  int vi = 0;
886  int comp_shift = 0;
887  for (i = 0; i < data->ne_component; i++) {
888  for (k = 0; k < data->ne_dof[i]; k++) {
889  /* Write for each block */
890  for (j = 0; j < num_blocks; j++) {
891  char vname[EXO_MAX_STR_LENGTH];
892  snprintf(vname, EXO_MAX_STR_LENGTH, "vals_elem_var%deb%d", vi + 1, j + 1);
893  int varid;
894  stat = nc_inq_varid(ncid, vname, &varid); NCERR(stat);
895 
896  int ne = blocks[j].num_elem;
897  double *vals = (double *)HECMW_malloc(sizeof(double) * ne);
898  int ei_idx;
899  for (ei_idx = 0; ei_idx < ne; ei_idx++) {
900  int ei = blocks[j].elem_indices[ei_idx];
901  vals[ei_idx] = data->elem_val_item[ei * data_tot_e + comp_shift + k];
902  }
903  size_t start[2] = {ts_idx, 0};
904  size_t count[2] = {1, (size_t)ne};
905  stat = nc_put_vara_double(ncid, varid, start, count, vals);
906  NCERR(stat);
907  HECMW_free(vals);
908  }
909  vi++;
910  }
911  comp_shift += data->ne_dof[i];
912  }
913  }
914 
915  /*========================================================================
916  * CLEANUP
917  *========================================================================*/
918  stat = nc_close(ncid);
919  NCERR(stat);
920 
921  /* Free element block index arrays */
922  for (j = 0; j < num_blocks; j++) {
923  HECMW_free(blocks[j].elem_indices);
924  }
925 
926  time_step_count++;
927 }
928 
929 /*============================================================================
930  * Public interface (same pattern as HECMW_vtk_output / HECMW_bin_vtk_output)
931  *============================================================================*/
932 
933 /* output_type=18 (EXODUS): All timesteps in one file */
934 void HECMW_exodus_output(struct hecmwST_local_mesh *mesh,
935  struct hecmwST_result_data *data,
936  char *outfile, char *outfile1,
937  HECMW_Comm VIS_COMM)
938 {
939  exodus_output(mesh, data, outfile, outfile1, VIS_COMM, 0);
940 }
941 
942 /* output_type=19 (STEP_EXODUS): One file per timestep */
943 void HECMW_exodus_step_output(struct hecmwST_local_mesh *mesh,
944  struct hecmwST_result_data *data,
945  char *outfile, char *outfile1,
946  HECMW_Comm VIS_COMM)
947 {
948  exodus_output(mesh, data, outfile, outfile1, VIS_COMM, 1);
949 }
950 
951 #else /* !WITH_NETCDF */
952 
953 #include <stdio.h>
954 #include "hecmw_struct.h"
955 #include "hecmw_result.h"
956 
957 /*============================================================================
958  * Stub functions when compiled without NetCDF support
959  *============================================================================*/
960 static void exodus_stub_warning(HECMW_Comm VIS_COMM)
961 {
962  int mynode;
963  HECMW_Comm_rank(VIS_COMM, &mynode);
964  if (mynode == 0) {
965  fprintf(stderr,
966  "WARNING: output_type=EXODUS requested but FrontISTR was not\n"
967  " compiled with NetCDF support.\n"
968  " Rebuild with -DWITH_NETCDF=ON (CMake) or\n"
969  " --with-netcdf (Makefile) to enable Exodus output.\n");
970  }
971 }
972 
973 void HECMW_exodus_output(struct hecmwST_local_mesh *mesh,
974  struct hecmwST_result_data *data,
975  char *outfile, char *outfile1,
976  HECMW_Comm VIS_COMM)
977 {
978  exodus_stub_warning(VIS_COMM);
979  (void)mesh; (void)data; (void)outfile; (void)outfile1;
980 }
981 
982 void HECMW_exodus_step_output(struct hecmwST_local_mesh *mesh,
983  struct hecmwST_result_data *data,
984  char *outfile, char *outfile1,
985  HECMW_Comm VIS_COMM)
986 {
987  exodus_stub_warning(VIS_COMM);
988  (void)mesh; (void)data; (void)outfile; (void)outfile1;
989 }
990 
991 #endif /* WITH_NETCDF */
if
if(!(yy_init))
Definition: hecmw_ablex.c:1823
HECMW_Comm_rank
int HECMW_Comm_rank(HECMW_Comm comm, int *rank)
Definition: hecmw_comm.c:18
HECMW_Comm_size
int HECMW_Comm_size(HECMW_Comm comm, int *size)
Definition: hecmw_comm.c:37
HECMW_FILENAME_LEN
#define HECMW_FILENAME_LEN
Definition: hecmw_config.h:72
HECMW_Comm
MPI_Comm HECMW_Comm
Definition: hecmw_config.h:30
HECMW_ctrl_make_subdir
int HECMW_ctrl_make_subdir(char *filename)
Definition: hecmw_control.c:2447
mesh
struct hecmwST_local_mesh * mesh
Definition: hecmw_repart.h:71
hecmw_etype.h
HECMW_exodus_output
void HECMW_exodus_output(struct hecmwST_local_mesh *mesh, struct hecmwST_result_data *data, char *outfile, char *outfile1, HECMW_Comm VIS_COMM)
Definition: hecmw_fstr_output_exo.c:973
HECMW_exodus_step_output
void HECMW_exodus_step_output(struct hecmwST_local_mesh *mesh, struct hecmwST_result_data *data, char *outfile, char *outfile1, HECMW_Comm VIS_COMM)
Definition: hecmw_fstr_output_exo.c:982
hecmw_fstr_output_exo.h
NULL
#define NULL
Definition: hecmw_io_nastran.c:30
hecmw_malloc.h
HECMW_free
#define HECMW_free(ptr)
Definition: hecmw_malloc.h:24
HECMW_malloc
#define HECMW_malloc(size)
Definition: hecmw_malloc.h:20
hecmw_result.h
hecmw_struct.h
hecmw_vis_comm_util.h
HECMW_vis_print_exit
void HECMW_vis_print_exit(char *var)
Definition: hecmw_vis_mem_util.c:21
hecmw_vis_mem_util.h
m_fstr::myrank
integer(kind=kint) myrank
PARALLEL EXECUTION.
Definition: m_fstr.f90:97
title
void title()
Definition: neu2fstr.cpp:52
data
CNFData data
Definition: neu_reporter.cpp:18
hecmwST_local_mesh
Definition: hecmw_struct.h:139
hecmwST_local_mesh::n_node
int n_node
Definition: hecmw_struct.h:161
hecmwST_local_mesh::node
double * node
Definition: hecmw_struct.h:170
hecmwST_local_mesh::elem_type
int * elem_type
Definition: hecmw_struct.h:191
hecmwST_local_mesh::elem_node_item
int * elem_node_item
Definition: hecmw_struct.h:196
hecmwST_local_mesh::global_node_ID
int * global_node_ID
Definition: hecmw_struct.h:168
hecmwST_local_mesh::elem_node_index
long long * elem_node_index
Definition: hecmw_struct.h:195
hecmwST_local_mesh::n_elem
int n_elem
Definition: hecmw_struct.h:184
hecmwST_local_mesh::global_elem_ID
int * global_elem_ID
Definition: hecmw_struct.h:190
hecmwST_result_data
Definition: hecmw_result.h:11