1 /*
2 Copyright (c) 1994 - 2010, Lawrence Livermore National Security, LLC.
3 LLNL-CODE-425250.
4 All rights reserved.
5
6 This file is part of Silo. For details, see silo.llnl.gov.
7
8 Redistribution and use in source and binary forms, with or without
9 modification, are permitted provided that the following conditions
10 are met:
11
12 * Redistributions of source code must retain the above copyright
13 notice, this list of conditions and the disclaimer below.
14 * Redistributions in binary form must reproduce the above copyright
15 notice, this list of conditions and the disclaimer (as noted
16 below) in the documentation and/or other materials provided with
17 the distribution.
18 * Neither the name of the LLNS/LLNL nor the names of its
19 contributors may be used to endorse or promote products derived
20 from this software without specific prior written permission.
21
22 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
23 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
24 LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
25 A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL LAWRENCE
26 LIVERMORE NATIONAL SECURITY, LLC, THE U.S. DEPARTMENT OF ENERGY OR
27 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
28 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
29 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
30 PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
31 LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
32 NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
33 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
34
35 This work was produced at Lawrence Livermore National Laboratory under
36 Contract No. DE-AC52-07NA27344 with the DOE. Neither the United
37 States Government nor Lawrence Livermore National Security, LLC nor
38 any of their employees, makes any warranty, express or implied, or
39 assumes any liability or responsibility for the accuracy,
40 completeness, or usefulness of any information, apparatus, product, or
41 process disclosed, or represents that its use would not infringe
42 privately-owned rights. Any reference herein to any specific
43 commercial products, process, or services by trade name, trademark,
44 manufacturer or otherwise does not necessarily constitute or imply its
45 endorsement, recommendation, or favoring by the United States
46 Government or Lawrence Livermore National Security, LLC. The views and
47 opinions of authors expressed herein do not necessarily state or
48 reflect those of the United States Government or Lawrence Livermore
49 National Security, LLC, and shall not be used for advertising or
50 product endorsement purposes.
51 */
52 #include <silo.h>
53 #include <stdio.h>
54 #include <stdlib.h>
55 #include <std.c>
56
57 void
add_edge(int nid1,int nid2,int * nedges,int * maxedges,int ** edges)58 add_edge(int nid1, int nid2, int *nedges, int *maxedges, int **edges)
59 {
60 if (*nedges == *maxedges)
61 {
62 *maxedges = (*maxedges) * 2 + 1;
63 *edges = (int *) realloc(*edges, *maxedges * 2 * sizeof(int));
64 }
65 (*edges)[2*(*nedges)+0] = nid1;
66 (*edges)[2*(*nedges)+1] = nid2;
67 *nedges = (*nedges) + 1;
68 }
69
70 int
have_edge(int nid1,int nid2,int nedges,const int * edges)71 have_edge(int nid1, int nid2, int nedges, const int *edges)
72 {
73 int i;
74 for (i = 0; i < nedges; i++)
75 {
76 if ((edges[2*i+0] == nid1 &&
77 edges[2*i+1] == nid2) ||
78 (edges[2*i+0] == nid2 &&
79 edges[2*i+1] == nid1))
80 return 1;
81 }
82 return 0;
83 }
84
85 #define HANDLE_EDGE(A,B) \
86 if (!have_edge(nodelist[nlidx+A],nodelist[nlidx+B],*nedges,*edges)) \
87 add_edge(nodelist[nlidx+A],nodelist[nlidx+B],nedges,&maxedges,edges)
88
89 void
build_edgelist(int nzones,int ndims,const int * nodelist,int lnodelist,int origin,int lo_off,int hi_off,const int * st,const int * sz,const int * sc,int nshapes,int * nedges,int ** edges)90 build_edgelist(int nzones, int ndims, const int *nodelist,
91 int lnodelist, int origin, int lo_off, int hi_off,
92 const int *st, const int *sz, const int *sc, int nshapes,
93 int *nedges, int **edges)
94 {
95 int nlidx = 0;
96 int maxedges = 0;
97 int shape;
98
99 *nedges = 0;
100 *edges = 0;
101 for (shape = 0; shape < nshapes; shape++)
102 {
103 int zonetype = st[shape];
104 int zonesize = sz[shape];
105 int zone;
106 for (zone = 0; zone < sc[shape]; zone++)
107 {
108 switch (zonetype)
109 {
110 case DB_ZONETYPE_HEX:
111 {
112 HANDLE_EDGE(0,1);
113 HANDLE_EDGE(0,3);
114 HANDLE_EDGE(0,4);
115 HANDLE_EDGE(1,2);
116 HANDLE_EDGE(1,5);
117 HANDLE_EDGE(2,3);
118 HANDLE_EDGE(2,6);
119 HANDLE_EDGE(3,7);
120 HANDLE_EDGE(4,5);
121 HANDLE_EDGE(4,7);
122 HANDLE_EDGE(5,6);
123 HANDLE_EDGE(6,7);
124 break;
125 }
126 }
127 nlidx += zonesize;
128 }
129 }
130 }
131
compare_nodes(const void * nap,const void * nbp)132 int compare_nodes(const void *nap, const void *nbp)
133 {
134 int na = *((int*)nap);
135 int nb = *((int*)nbp);
136 if (na < nb) return -1;
137 else if (na > nb) return 1;
138 return 0;
139 }
140
141 void
add_face(int nid1,int nid2,int nid3,int nid4,int * nfaces,int * maxfaces,int ** faces)142 add_face(int nid1, int nid2, int nid3, int nid4, int *nfaces, int *maxfaces, int **faces)
143 {
144 int i;
145 int fnid[4];
146 fnid[0] = nid1;
147 fnid[1] = nid2;
148 fnid[2] = nid3;
149 fnid[3] = nid4;
150 qsort(fnid, sizeof(int), 4, compare_nodes);
151 if (*nfaces == *maxfaces)
152 {
153 *maxfaces = (*maxfaces) * 2 + 1;
154 *faces = (int *) realloc(*faces, *maxfaces * 4 * sizeof(int));
155 }
156 for (i = 0; i < 4; i++)
157 (*faces)[4*(*nfaces)+i] = fnid[i];
158 *nfaces = (*nfaces) + 1;
159 }
160
161 int
have_face(int nid1,int nid2,int nid3,int nid4,int nfaces,const int * faces)162 have_face(int nid1, int nid2, int nid3, int nid4, int nfaces, const int *faces)
163 {
164 int i,j;
165 int fnid[4];
166 fnid[0] = nid1;
167 fnid[1] = nid2;
168 fnid[2] = nid3;
169 fnid[3] = nid4;
170 qsort(fnid, sizeof(int), 4, compare_nodes);
171 for (i = 0; i < nfaces; i++)
172 {
173 int allmatch = 1;
174 for (j = 0; j < 4; j++)
175 {
176 if (faces[4*i+j] != fnid[j])
177 {
178 allmatch = 0;
179 break;
180 }
181 }
182 if (allmatch) return 1;
183 }
184 return 0;
185 }
186
187 #define HANDLE_FACE(A,B,C,D) \
188 if (!have_face(nodelist[nlidx+A],nodelist[nlidx+B], \
189 nodelist[nlidx+C],nodelist[nlidx+D],*nfaces,*faces)) \
190 add_face(nodelist[nlidx+A],nodelist[nlidx+B], \
191 nodelist[nlidx+C],nodelist[nlidx+D],nfaces,&maxfaces,faces)
192
193 void
build_facelist(int nzones,int ndims,const int * nodelist,int lnodelist,int origin,int lo_off,int hi_off,const int * st,const int * sz,const int * sc,int nshapes,int * nfaces,int ** faces)194 build_facelist(int nzones, int ndims, const int *nodelist,
195 int lnodelist, int origin, int lo_off, int hi_off,
196 const int *st, const int *sz, const int *sc, int nshapes,
197 int *nfaces, int **faces)
198 {
199 int nlidx = 0;
200 int maxfaces = 0;
201 int shape;
202
203 *nfaces = 0;
204 *faces = 0;
205 for (shape = 0; shape < nshapes; shape++)
206 {
207 int zonetype = st[shape];
208 int zonesize = sz[shape];
209 int zone;
210 for (zone = 0; zone < sc[shape]; zone++)
211 {
212 switch (zonetype)
213 {
214 case DB_ZONETYPE_HEX:
215 {
216 HANDLE_FACE(0,1,5,4);
217 HANDLE_FACE(0,3,2,1);
218 HANDLE_FACE(0,4,7,3);
219 HANDLE_FACE(1,2,6,5);
220 HANDLE_FACE(2,3,7,6);
221 HANDLE_FACE(4,5,6,7);
222 break;
223 }
224 }
225 nlidx += zonesize;
226 }
227 }
228 }
229
230 int
main(int argc,char * argv[])231 main(int argc, char *argv[])
232 {
233 DBfile *dbfile = NULL;
234 char *coordnames[3];
235 float *coords[3];
236 float x[64], y[64], z[64];
237 int shapesize[1];
238 int shapecnt[1];
239 DBfacelist *facelist = NULL;
240 int matnos[1], matlist[1], dims[3];
241 int i, j, k, len;
242 float evar2d[2*16], evar3d[3*64], fvar3d[3*64];
243 int driver = DB_PDB;
244 char *filename = "efcentering.silo";
245 int layer, zone;
246 int nodelist2[9*4] = {0,1,5,4, 1,2,6,5, 2,3,7,6,
247 4,5,9,8, 5,6,10,9, 6,7,11,10,
248 8,9,13,12, 9,10,14,13, 10,11,15,14};
249 int st2 = DB_ZONETYPE_QUAD;
250 int ss2 = 4;
251 int sc2 = 9;
252 int nodelist3[27*8];
253 int st3 = DB_ZONETYPE_HEX;
254 int ss3 = 8;
255 int sc3 = 27;
256
257 int nedges;
258 int *edges;
259 int nfaces;
260 int *faces;
261 int ndims;
262 int show_all_errors = FALSE;
263
264 /* Parse command-line */
265 for (i=1; i<argc; i++) {
266 if (!strncmp(argv[i], "DB_PDB",6)) {
267 driver = StringToDriver(argv[i]);
268 } else if (!strncmp(argv[i], "DB_HDF5", 7)) {
269 driver = StringToDriver(argv[i]);
270 } else if (!strcmp(argv[i], "show-all-errors")) {
271 show_all_errors = 1;
272 } else if (argv[i][0] != '\0') {
273 fprintf(stderr, "%s: ignored argument `%s'\n", argv[0], argv[i]);
274 }
275 }
276
277 DBShowErrors(show_all_errors?DB_ALL_AND_DRVR:DB_ABORT, NULL);
278 dbfile = DBCreate(filename, DB_CLOBBER, DB_LOCAL, "edge and face centered data", driver);
279
280 coordnames[0] = "xcoords";
281 coordnames[1] = "ycoords";
282 coordnames[2] = "zcoords";
283
284 dims[0] = 4;
285 dims[1] = 4;
286 dims[2] = 4;
287 for (k = 0; k < 4; k++)
288 {
289 for (j = 0; j < 4; j++)
290 {
291 for (i = 0; i < 4; i++)
292 {
293 x[k*4*4+j*4+i] = (float) i;
294 y[k*4*4+j*4+i] = (float) j;
295 z[k*4*4+j*4+i] = (float) k;
296 evar2d[0*16+j*4+i] = (float) i;
297 evar2d[1*16+j*4+i] = (float) j;
298 evar3d[0*64+k*4*4+j*4+i] = (float) i;
299 evar3d[1*64+k*4*4+j*4+i] = (float) j;
300 evar3d[2*64+k*4*4+j*4+i] = (float) k;
301 fvar3d[0*64+k*4*4+j*4+i] = (float) 10*i;
302 fvar3d[1*64+k*4*4+j*4+i] = (float) 100*j;
303 fvar3d[2*64+k*4*4+j*4+i] = (float) 1000*k;
304 }
305 }
306 }
307
308 coords[0] = x;
309 coords[1] = y;
310 coords[2] = z;
311
312 /* build 3d zonelist by layering 2d zonelist */
313 for (layer = 0; layer < 3; layer++)
314 {
315 for (zone = 0; zone < 9; zone++)
316 {
317 nodelist3[layer*9*8+zone*8+0] = nodelist2[zone*4+0]+(layer+1)*16;
318 nodelist3[layer*9*8+zone*8+1] = nodelist2[zone*4+0]+layer*16;
319 nodelist3[layer*9*8+zone*8+2] = nodelist2[zone*4+1]+layer*16;
320 nodelist3[layer*9*8+zone*8+3] = nodelist2[zone*4+1]+(layer+1)*16;
321 nodelist3[layer*9*8+zone*8+4] = nodelist2[zone*4+3]+(layer+1)*16;
322 nodelist3[layer*9*8+zone*8+5] = nodelist2[zone*4+3]+layer*16;
323 nodelist3[layer*9*8+zone*8+6] = nodelist2[zone*4+2]+layer*16;
324 nodelist3[layer*9*8+zone*8+7] = nodelist2[zone*4+2]+(layer+1)*16;
325 }
326 }
327
328 DBPutQuadmesh(dbfile, "qmesh2", (DBCAS_t) coordnames, coords, dims, 2, DB_FLOAT, DB_NONCOLLINEAR, 0);
329 DBPutQuadmesh(dbfile, "qmesh3", (DBCAS_t) coordnames, coords, dims, 3, DB_FLOAT, DB_NONCOLLINEAR, 0);
330 DBPutQuadvar1(dbfile, "qevar2", "qmesh2", evar2d, dims, 2, 0, 0, DB_FLOAT, DB_EDGECENT, 0);
331 DBPutQuadvar1(dbfile, "qevar3", "qmesh3", evar3d, dims, 3, 0, 0, DB_FLOAT, DB_EDGECENT, 0);
332 DBPutQuadvar1(dbfile, "qfvar3", "qmesh3", fvar3d, dims, 3, 0, 0, DB_FLOAT, DB_FACECENT, 0);
333
334 DBPutUcdmesh(dbfile, "umesh2", 2, (DBCAS_t) coordnames, coords, 16, 9, "um2zl", 0, DB_FLOAT, 0);
335 DBPutUcdmesh(dbfile, "umesh3", 3, (DBCAS_t) coordnames, coords, 64, 27, "um3zl", 0, DB_FLOAT, 0);
336 DBPutZonelist2(dbfile, "um2zl", 9, 2, nodelist2, ss2*sc2, 0, 0, 0, &st2, &ss2, &sc2, 1, 0);
337 DBPutZonelist2(dbfile, "um3zl", 27, 3, nodelist3, ss3*sc3, 0, 0, 0, &st3, &ss3, &sc3, 1, 0);
338
339 /* Only reason we build an edgelist is so we know the number of unique edges in the mesh */
340 build_edgelist(27, 3, nodelist3, ss3*sc3, 0, 0, 0, &st3, &ss3, &sc3, 1, &nedges, &edges);
341 for (i = 0; i < nedges; i++)
342 evar3d[i] = i;
343 DBPutUcdvar1(dbfile, "uevar3", "umesh3", evar3d, nedges, 0, 0, DB_FLOAT, DB_EDGECENT, 0);
344 ndims = 2;
345 dims[0] = nedges;
346 dims[1] = 2;
347 DBWrite(dbfile, "edges", edges, dims, ndims, DB_INT);
348 free(edges);
349
350 /* Only reason we build a facelist is so we know the number of unique faces in the mesh */
351 build_facelist(27, 3, nodelist3, ss3*sc3, 0, 0, 0, &st3, &ss3, &sc3, 1, &nfaces, &faces);
352 for (i = 0; i < nfaces; i++)
353 fvar3d[i] = i;
354 DBPutUcdvar1(dbfile, "ufvar3", "umesh3", fvar3d, nfaces, 0, 0, DB_FLOAT, DB_FACECENT, 0);
355 dims[0] = nfaces;
356 dims[1] = 4;
357 DBWrite(dbfile, "faces", faces, dims, ndims, DB_INT);
358 free(faces);
359
360 DBClose(dbfile);
361
362 CleanupDriverStuff();
363 return (0);
364 }
365