1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <math.h>
5 #ifdef _WIN32
6 # include <io.h>
7 # define unlink _unlink
8 #else
9 # include <unistd.h>
10 #endif
11 #include "cgnslib.h"
12
13 #ifndef CGNSTYPES_H
14 # define cgsize_t int
15 #endif
16 #ifndef CGNS_ENUMT
17 # define CGNS_ENUMT(e) e
18 # define CGNS_ENUMV(e) e
19 #endif
20
21 int CellDim = 3, PhyDim = 3;
22 int cgfile, cgbase, cgzone, cggrid, cgsect, cgsol, cgfld;
23
24 cgsize_t size[9] = {5, 5, 5, 4, 4, 4, 0, 0, 0};
25 int rind[6] = {0, 0, 0, 0, 1, 1};
26
27 #define NUM_I (size[0] + rind[0] + rind[1])
28 #define NUM_J (size[1] + rind[2] + rind[3])
29 #define NUM_K (size[2] + rind[4] + rind[5])
30
31 #define INDEX(I,J,K) (int)((I) + NUM_I * ((J) + NUM_J * (K)))
32
33 cgsize_t num_coord;
34 float *xcoord, *ycoord, *zcoord;
35 cgsize_t *nmap;
36
37 cgsize_t num_element, *elements;
38
39 float *solution;
40
compute_coord(int n,int i,int j,int k)41 static void compute_coord (int n, int i, int j, int k)
42 {
43 xcoord[n] = (float)(i - rind[0]);
44 ycoord[n] = (float)(j - rind[2]);
45 zcoord[n] = (float)(k - rind[4]);
46 }
47
compute_element(int ne,int i,int j,int k)48 static void compute_element (int ne, int i, int j, int k)
49 {
50 int n = ne << 3;
51
52 elements[n++] = nmap[INDEX(i, j, k)];
53 elements[n++] = nmap[INDEX(i+1, j, k)];
54 elements[n++] = nmap[INDEX(i+1, j+1, k)];
55 elements[n++] = nmap[INDEX(i, j+1, k)];
56 elements[n++] = nmap[INDEX(i, j, k+1)];
57 elements[n++] = nmap[INDEX(i+1, j, k+1)];
58 elements[n++] = nmap[INDEX(i+1, j+1, k+1)];
59 elements[n] = nmap[INDEX(i, j+1, k+1)];
60 }
61
main(int argc,char * argv[])62 int main (int argc, char *argv[])
63 {
64 int n, nn, i, j, k, ni, nj, nk;
65 cgsize_t dims[3];
66 int grind[2], erind[2];
67
68 if (argc > 1) {
69 n = 0;
70 if (argv[1][n] == '-') n++;
71 if (argv[1][n] == 'a' || argv[1][n] == 'A')
72 nn = CG_FILE_ADF;
73 else if (argv[1][n] == 'h' || argv[1][n] == 'H')
74 nn = CG_FILE_HDF5;
75 else {
76 fprintf(stderr, "unknown option\n");
77 exit (1);
78 }
79 if (cg_set_file_type(nn))
80 cg_error_exit();
81 }
82
83 num_coord = NUM_I * NUM_J * NUM_K;
84 num_element = (NUM_I - 1) * (NUM_J - 1) * (NUM_K - 1);
85 xcoord = (float *) malloc((size_t)(4 * num_coord * sizeof(float)));
86 nmap = (cgsize_t *) malloc((size_t)((num_coord + 8 * num_element) * sizeof(cgsize_t)));
87 if (NULL == xcoord || NULL == nmap) {
88 fprintf(stderr, "malloc failed for data\n");
89 exit(1);
90 }
91 ycoord = xcoord + num_coord;
92 zcoord = ycoord + num_coord;
93 solution = zcoord + num_coord;
94 elements = nmap + num_coord;
95
96 for (n = 0; n < num_coord; n++)
97 solution[n] = (float)n;
98
99 unlink("rind.cgns");
100 if (cg_open("rind.cgns", CG_MODE_WRITE, &cgfile))
101 cg_error_exit();
102
103 /*--- structured grid with rind ---*/
104
105 printf ("writing structured base with rind\n");
106 fflush (stdout);
107
108 for (n = 0, k = 0; k < NUM_K; k++) {
109 for (j = 0; j < NUM_J; j++) {
110 for (i = 0; i < NUM_I; i++) {
111 compute_coord(n++, i, j, k);
112 }
113 }
114 }
115
116 if (cg_base_write(cgfile, "Structured", CellDim, PhyDim, &cgbase) ||
117 cg_goto(cgfile, cgbase, "end") ||
118 cg_dataclass_write(CGNS_ENUMV( NormalizedByUnknownDimensional )) ||
119 cg_zone_write(cgfile, cgbase, "Zone", size, CGNS_ENUMV( Structured ), &cgzone))
120 cg_error_exit();
121
122 /* can't use cg_coord_write to write rind coordinates
123 need to use cg_grid_write to create the node, cg_goto to set
124 position at the node, then write rind and coordinates as an array */
125
126 dims[0] = NUM_I;
127 dims[1] = NUM_J;
128 dims[2] = NUM_K;
129
130 if (cg_grid_write(cgfile, cgbase, cgzone, "GridCoordinates", &cggrid) ||
131 cg_goto(cgfile, cgbase, "Zone_t", cgzone,
132 "GridCoordinates_t", cggrid, "end") ||
133 cg_rind_write(rind) ||
134 cg_array_write("CoordinateX", CGNS_ENUMV( RealSingle ), 3, dims, xcoord) ||
135 cg_array_write("CoordinateY", CGNS_ENUMV( RealSingle ), 3, dims, ycoord) ||
136 cg_array_write("CoordinateZ", CGNS_ENUMV( RealSingle ), 3, dims, zcoord))
137 cg_error_exit();
138
139 /* a similar technique is used for the solution with rind,
140 but we use cg_field_write instead of cg_array_write
141 and the solution dimensions come from the zone sizes */
142
143 if (cg_sol_write(cgfile, cgbase, cgzone, "VertexSolution",
144 CGNS_ENUMV( Vertex ), &cgsol) ||
145 cg_goto(cgfile, cgbase, "Zone_t", cgzone,
146 "FlowSolution_t", cgsol, "end") ||
147 cg_rind_write(rind) ||
148 cg_field_write(cgfile, cgbase, cgzone, cgsol, CGNS_ENUMV( RealSingle ),
149 "Density", solution, &cgfld))
150 cg_error_exit();
151
152 /*--- unstructured with rind ---*/
153
154 printf ("writing unstructured base with rind\n");
155 fflush (stdout);
156
157 /* rind here has dimension rind[2], so need to put all the
158 rind coordinates at the beginning and/or end of the array.
159 Just for grins, I'll put some at both ends, although
160 it's probably best to put the rind coordinates at the end.
161 I'll use the nmap array for building elements */
162
163 ni = (int)(size[0] + rind[0]);
164 nj = (int)(size[1] + rind[2]);
165 nk = (int)(size[2] + rind[4]);
166
167 for (n = 0, i = 0; i < rind[0]; i++) {
168 for (k = 0; k < NUM_K; k++) {
169 for (j = 0; j < NUM_J; j++) {
170 compute_coord (n, i, j, k);
171 nn = INDEX(i, j, k);
172 nmap[nn] = ++n;
173 }
174 }
175 }
176 for (j = 0; j < rind[2]; j++) {
177 for (k = 0; k < NUM_K; k++) {
178 for (i = rind[0]; i < ni; i++) {
179 compute_coord (n, i, j, k);
180 nn = INDEX(i, j, k);
181 nmap[nn] = ++n;
182 }
183 }
184 }
185 for (k = 0; k < rind[4]; k++) {
186 for (j = rind[2]; j < nj; j++) {
187 for (i = rind[0]; i < ni; i++) {
188 compute_coord (n, i, j, k);
189 nn = INDEX(i, j, k);
190 nmap[nn] = ++n;
191 }
192 }
193 }
194 grind[0] = n;
195
196 for (k = rind[4]; k < nk; k++) {
197 for (j = rind[2]; j < nj; j++) {
198 for (i = rind[0]; i < ni; i++) {
199 compute_coord (n, i, j, k);
200 nn = INDEX(i, j, k);
201 nmap[nn] = ++n;
202 }
203 }
204 }
205 grind[1] = (int)(num_coord - n);
206
207 for (i = ni; i < NUM_I; i++) {
208 for (k = 0; k < NUM_K; k++) {
209 for (j = 0; j < NUM_J; j++) {
210 compute_coord (n, i, j, k);
211 nn = INDEX(i, j, k);
212 nmap[nn] = ++n;
213 }
214 }
215 }
216 for (j = nj; j < NUM_J; j++) {
217 for (k = 0; k < NUM_K; k++) {
218 for (i = rind[0]; i < ni; i++) {
219 compute_coord (n, i, j, k);
220 nn = INDEX(i, j, k);
221 nmap[nn] = ++n;
222 }
223 }
224 }
225 for (k = nk; k < NUM_K; k++) {
226 for (j = rind[2]; j < nj; j++) {
227 for (i = rind[0]; i < ni; i++) {
228 compute_coord (n, i, j, k);
229 nn = INDEX(i, j, k);
230 nmap[nn] = ++n;
231 }
232 }
233 }
234
235 /* rind elements are like the coordinates, they need to go
236 at the beginning and/or end of the element array, although
237 at the end is probably best */
238
239 for (n = 0, i = 0; i < rind[0]; i++) {
240 for (k = 0; k < NUM_K - 1; k++) {
241 for (j = 0; j < NUM_J - 1; j++) {
242 compute_element(n++, i, j, k);
243 }
244 }
245 }
246 for (j = 0; j < rind[2]; j++) {
247 for (k = 0; k < NUM_K - 1; k++) {
248 for (i = rind[0]; i < ni - 1; i++) {
249 compute_element(n++, i, j, k);
250 }
251 }
252 }
253 for (k = 0; k < rind[4]; k++) {
254 for (j = rind[2]; j < nj - 1; j++) {
255 for (i = rind[0]; i < ni - 1; i++) {
256 compute_element(n++, i, j, k);
257 }
258 }
259 }
260 erind[0] = n;
261
262 for (k = rind[4]; k < nk - 1; k++) {
263 for (j = rind[2]; j < nj - 1; j++) {
264 for (i = rind[0]; i < ni - 1; i++) {
265 compute_element(n++, i, j, k);
266 }
267 }
268 }
269 erind[1] = (int)(num_element - n);
270
271 for (i = ni - 1; i < NUM_I - 1; i++) {
272 for (k = 0; k < NUM_K - 1; k++) {
273 for (j = 0; j < NUM_J - 1; j++) {
274 compute_element(n++, i, j, k);
275 }
276 }
277 }
278 for (j = nj - 1; j < NUM_J - 1; j++) {
279 for (k = 0; k < NUM_K - 1; k++) {
280 for (i = rind[0]; i < ni - 1; i++) {
281 compute_element(n++, i, j, k);
282 }
283 }
284 }
285 for (k = nk - 1; k < NUM_K - 1; k++) {
286 for (j = rind[2]; j < nj - 1; j++) {
287 for (i = rind[0]; i < ni - 1; i++) {
288 compute_element(n++, i, j, k);
289 }
290 }
291 }
292
293 /* create base, zone and write coordinates.
294 As for the structured case, the rind coordinates
295 and elements are not included in the zone totals */
296
297 dims[0] = num_coord - grind[0] - grind[1];
298 dims[1] = num_element - erind[0] - erind[1];
299 dims[2] = 0;
300
301 if (cg_base_write(cgfile, "Unstructured", CellDim, PhyDim, &cgbase) ||
302 cg_goto(cgfile, cgbase, "end") ||
303 cg_dataclass_write(CGNS_ENUMV( NormalizedByUnknownDimensional )) ||
304 cg_zone_write(cgfile, cgbase, "Zone", dims, CGNS_ENUMV( Unstructured ), &cgzone) ||
305 cg_grid_write(cgfile, cgbase, cgzone, "GridCoordinates", &cggrid) ||
306 cg_goto(cgfile, cgbase, "Zone_t", cgzone,
307 "GridCoordinates_t", cggrid, "end") ||
308 cg_rind_write(grind) ||
309 cg_array_write("CoordinateX", CGNS_ENUMV( RealSingle ), 1, &num_coord, xcoord) ||
310 cg_array_write("CoordinateY", CGNS_ENUMV( RealSingle ), 1, &num_coord, ycoord) ||
311 cg_array_write("CoordinateZ", CGNS_ENUMV( RealSingle ), 1, &num_coord, zcoord))
312 cg_error_exit();
313
314 /* to write the elements with rind elements,
315 write all the elements, then use goto to add rind */
316
317 if (cg_section_write(cgfile, cgbase, cgzone, "Elements", CGNS_ENUMV( HEXA_8 ),
318 1, num_element, 0, elements, &cgsect) ||
319 cg_goto(cgfile, cgbase, "Zone_t", cgzone,
320 "Elements_t", cgsect, "end") ||
321 cg_rind_write(erind))
322 cg_error_exit();
323
324 /* write solution a vertices with rind */
325
326 if (cg_sol_write(cgfile, cgbase, cgzone, "VertexSolution",
327 CGNS_ENUMV( Vertex ), &cgsol) ||
328 cg_goto(cgfile, cgbase, "Zone_t", cgzone,
329 "FlowSolution_t", cgsol, "end") ||
330 cg_rind_write(grind) ||
331 cg_field_write(cgfile, cgbase, cgzone, cgsol, CGNS_ENUMV( RealSingle ),
332 "Density", solution, &cgfld))
333 cg_error_exit();
334
335 cg_close(cgfile);
336 return 0;
337 }
338
339